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INRTODCUCTION 


This  postdoctoral  traineeship  grant  (DAMD17-03-1-0657,  entitled  “Multiple  Aperture  Radiation 
Therapy  for  Breast  Cancer”)  was  awarded  for  the  period  of  Nov.  1,  2003  to  Oct.  31,  2005,  and 
was  transferred  to  the  current  principal  investigator  (PI)  on  Nov.  1,  2004  and  extended  without 
funds  to  November  2006.  The  goal  of  this  project  is  to  develop  a  novel  radiation  treatment 
technique,  multiple  aperture  radiation  therapy  (MART),  as  a  candidate  modality  for  treating 
breast  cancer.  Under  the  generous  support  from  the  U.S.  Army  Medical  Research  and  Materiel 
Command  (USAMRMC),  the  PI  has  gained  a  tremendous  amount  of  knowledge  on  breast  cancer 
and  breast  cancer  management.  The  support  has  also  made  it  possible  for  the  PI  to  contribute 
significantly  to  breast  cancer  research.  A  number  of  conference  and  refereed  publications  have 
been  resulted  from  the  support.  Two  annual  reports  have  been  previously  submitted  and 
reviewed  positively,  and  this  one  is  the  final  comprehensive  report  to  summarize  all  the  scientific 
results  and  achievements  made  during  the  work  of  the  project. 


BODY 


A.  Background 

Breast  cancer  is  the  most  frequently  diagnosed  cancer  and  the  second  leading  cause  of 
cancer  death  in  women.  The  American  Cancer  Society  estimates  a  total  of  212,920  new  female 
breast  cancers  and  an  estimated  40,970  deaths  from  them  in  the  US  during  2006  [1],  For  patients 
with  early  breast  cancer,  radiation  therapy  is  commonly  used  in  the  treatment.  Conventional 
radiotherapy  for  breast  cancer  utilizes  two  opposed  tangential  fields  (OTF)  with  either  uniform 
or  wedged  photon  beams  [2-5].  While  being  effective  the  conventional  method  presents 
problems  related  to  breast  dose  inhomogeneity  and  relatively  high  doses  to  the  ipsilateral  lung 
and  heart.  Modem  therapeutic  advance  in  Intensity  Modulated  Radiation  Therapy  (IMRT) 
provides  unprecedented  means  to  deliver  3D-dose  distributions  with  superior  tumor  conformality 
and  normal  tissue  spare,  and  can  potentially  overcome  these  problems  9'14.  However,  there  are 
several  practical  issues  in  breast  IMRT  treatment:  the  planning,  delivery  and  quality  assurance 
(QA)  processes  add  significant  overhead  as  compared  to  the  standard  OTF  and  make  it  clinically 
less  useful.  For  instance,  current  IMRT  treatment  planning  requires  physicians  to  segment 
explicitly  the  target  volume,  which  takes  in  average  about  20  minutes  per  patient  and  increases 
the  cost  of  health  care.  The  procedure  of  dosimetric  verification  of  a  complicate  shape  intensity- 
modulated  beam  is  yet  not  well  established  and  needs  great  effort  of  quality  control.  In  addition, 
the  IMRT  optimization  algorithm  is  not  dedicated  to  breast  cancer  and  may  not  be  the  most 
efficient  scheme.  In  short,  an  effective  method  to  plan  and  deliver  IMRT  breast  treatment  is 
highly  desired  in  order  for  tens  of  thousands  of  breast  cancer  patients  to  benefit  from  the  state-of- 
the-art  technology.  For  this  purpose,  we  have  proposed  the  MART  approach  in  this  USAMRMC- 
supported  project  to  combine  the  planning  simplicity  of  OFT  with  the  precise  dose  control  of 
IMRT.  To  achieve  this,  MART  has  been  designed  by  adding  finite  beam  segments  to  a  properly 
selected  tangential  field,  so  that  the  dose  inhomogeneity  resulted  from  OTF  is  efficiently 
removed  and  the  total  dose  are  better  confined  to  the  target  volume.  During  the  development  and 
implementation  of  the  proposed  MART  technique,  we  also  devoted  ourselves  into  solving  a 
practical  issue  in  breast  treatment  related  to  the  breathing  motion.  The  breathing  motion  affects 
the  position  of  the  breast,  lung  and  heart  during  the  treatment,  and  thus  introduces  extra 
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geometric  uncertainties.  Conventionally,  these  uncertainties  have  to  be  accommodated  by 
enlarging  the  treatment  margin,  resulting  in  relatively  high  dose  to  the  normal  tissues.  To  reduce 
the  margin  and  improve  the  efficacy  of  MART,  we  have  extensively  studied  a  variety  of 
techniques  related  to  the  motion  compensation  or  four-dimensional  (4D)  treatment.  The  details 
are  summarized  below. 

B.  MART  Development  (Nov  2003  ~  Nov  2004) 

B.l.  Manual  MART  planning 

We  have  developed  and 
clinically  implemented  a  manual 
MART  planning  procedure  for 
breast  irradiation  at  Stanford 
University  Hospital  in  2004.  The 
manual  treatment  planning  was 
done  with  a  3D  treatment  planning 
system  (FOCUS™,  Computerized 
Medical  System,  St.  Louis,  MO).  It 
started  with  a  standard  OTF  plan. 

Basically,  patients  underwent 
computed  tomography  (CT)  in  the 
conventional  treatment  position 
supported  by  an  Alpha  Cradle 
immobilization  device  (Smithers 
Medical  Products,  Tallmadge, 

Ohio).  Radiopaque  markers  were 
placed  on  the  patients’  chest  to 
indicate  the  medial  and  lateral 
borders  of  the  palpable  breast  tissue 
and  the  location  of  the  lumpectomy 
scar.  The  radiation-sensitive 
structures  included  the  left  and  right 
lungs,  the  heart,  and  the 
contralateral  breast.  For  the  purpose 
of  quantitative  study,  the  contours  of  the  skin,  target  volume  and  the  sensitive  structures  were 
outlined  using  the  segmentation  tools  provided  by  the  virtual  simulation  workstation 
(AcQSim™,  Philips  Medical  System,  Cleveland,  OH).  It  is,  however,  not  required  to  outline  the 
structures  in  general  MART  treatment.  The  tangential  fields  were  determined  by  the  routine 
virtual  simulation  procedure  performed  on  an  AcQSim  workstation.  The  fields  may  be  adjusted 
at  the  stage  of  treatment  planning  according  to  the  actual  treatment  objective  for  each  patient 
with  considerations  concerning  tumor  bed  coverage,  in-field  lung  and  cardiac  volume,  if  left 
breast  irradiation.  Figure  la  shows  an  example  of  the  OTF  setup  for  the  treatment  of  a  left  breast 
cancer  patient.  After  this,  we  proceeded  to  introduce  an  additional  MLC  field  segment  to  one  or 
both  beam  directions  to  boost  the  “cold”  region(s)  under  the  guidance  of  dose  distributions  in  the 
plane  perpendicular  to  the  incident  beam  direction.  Figures  1  (d)-(i)  show  the  three  segments  of 
the  lateral  and  the  medial  fields  for  a  MART  treatment.  The  weights  and  MLC  apertures  of  the 
segments  were  adjusted  manually  to  achieve  a  uniform  dose  distribution.  Our  experience 
indicated  that,  for  intermediately  complex  cases,  it  was  often  sufficient  to  introduce  one  or  two 


Figure  1  Standard  tangential  field  arrangement  for  treatment  of  a 
left  breast  cancer  patient  (top  row).  The  middle  and  bottom  rows  are 
the  MLC  shapes  of  the  three  segments  of  the  medial  and  lateral 
MART  fields,  respectively.  A  physical  wedge  of  30°  was  used  in 
the  lateral  field. 
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additional  segments  to  the  original  opposed  tangential  fields.  For  complex  cases,  two  or  three 
additional  segments  were  frequently  used.  A  segment  with  ipsilateral  lung  or  heart  blocked  by 
MLC  (the  third  segment  in  either  medial  or  lateral  MART  field)  was  also  helpful  in  reducing  the 
dose  to  these  structures. 

B.2.  Optimization  for  MART 

We  have  developed  an  effective  aperture -based  optimization  algorithm  for  breast 
irradiation:  a  genetic  algorithm  (GA)  for  breast  MART  inverse  planning.  It  optimizes 
simultaneously  the  shapes  of  a  pre-specified  number  of  segments  and  their  corresponding 
weights.  The  GA  encodes  potential  solutions  in  chromosome-like  structures  and  applies  different 
recombination  operators  (crossover  and  mutation)  to  explore  the  search  space.  The  system 
variables,  which  include  the  weights  of  the  segments  and  the  leaf  positions  defining  the  shapes  of 
the  segments  are  encoded  in  three  chromosome-like  structures  for  each  potential  solution:  the  1st 
one  for  the  weights  of  the  segments,  the  2nd  one  for  the  positions  of  left  bank  of  MLC  leaves, 
and  3rd  one  for  the  right  bank.  All  encodings  are  realized  in  integer  format.  The  quality  of  a 
solution  is  evaluated  according  to  its  fitness,  which  is  defined  as  the  inverse  of  the  usual 
quadratic  objective  value.  The  chance  of  a  particular  solution  to  be  part  of  the  next  generation 
population  is  proportional  to  its  survival  probability,  defined  as  the  ratio  of  its  fitness  and  the 
total  fitness  of  the  population.  A  more  fit  solution  has  a  higher  probability  of  being  selected  into 
the  next  generation.  The  new  population  is  then  selected  by  simulating  the  spinning  of  a  suitable 
roulette  wheel,  N  times,  where  N  equals  the  number  of  solutions  in  the  current  generation.  The 
selection  process  is  followed  by  crossover  and  mutation  operations,  where  the  potential  solutions 
interchange  information  that  usually  leads  to  improved  individuals.  In  addition  to  the  fitness- 
based  selection  process,  we  allow  the  best  member  of  the  current  population  to  be  automatically 
copied  into  the  next  generation.  This  process  called  elitism  leads  to  a  faster  convergence  and 
keeps  track  of  the  best  solutions  obtained  in  each  iteration.  The  segment-based  optimized  plans 
improved  significantly  the  target  dose  uniformity  in  comparison  with  the  standard  OTF  plans. 
The  overall  planning  and  treatment  delivery  overhead  of  the  approach  is  significantly  reduced 
compared  to  the  conventional  beamlet-based  IMRT  breast  irradiation.  For  all  cases  (5  left  and  5 
right  breast  cancer  patients)  we  have  tested,  it  was  found  that  3~10  segments  per  tangential  beam 
were  sufficient  to  ensure  highly  homogeneous  doses  within  the  target.  Our  results  also  revealed 
that  the  maximum  target  dose  could  easily  be  reduced  from  109%~1 17%  in  conventional  OTF  to 
106%  to  112%.  The  volume  receiving  high  dose  irradiation  in  the  breast  target  was  also 
markedly  reduced.  It  was  also  possible  to  use  segment-based  optimization  to  reduce  the  dose  to 
the  ipsilateral  lung/heart. 

Another  significant  technique  that  we  have  developed  for  the  MART  inverse  planning  is 
to  purposely  modulate  the  penalty  on  individual  voxel  level  based  on  the  a  priori  dosimetric 
information  of  the  dose  optimization  system.  This  is  based  on  the  fact  that  the  voxels  within  a 
structure  are  not  identical  in  complying  with  their  dosiemtric  goals  and  there  exists  strong  intra- 
structural  competition  in  these  voxels.  Inverse  planning  objective  function  should  not  only 
balance  the  competing  objectives  of  different  structures  but  also  that  of  the  individual  voxels 
within  various  structures.  We  quantify  the  degree  for  a  voxel  to  achieve  its  dosimetric  goal  by 
introducing  the  concept  of  dosimetric  capability  for  each  voxel  in  a  target  or  sensitive  structure. 
To  the  contrary,  conventional  inverse  planning  algorithms  treat  all  voxels  within  a  target  or 
sensitive  structure  equally  and  use  structure  specific  prescriptions  and  weighting  factors  as 
system  parameters. 
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C.  Impact  of  MART  (Nov  2004  ~  Nov  2006) 


For  the  evaluation  of  the  MART  technique,  two  types  of  plans  were  generated  and 
compared  for  35  patient  enrolled  in  the  study.  One  was  the  standard  OTF  plan  and  the  other  one 
was  the  MART  plan.  Standard  plan  involved  a  medial  and  lateral  tangential  field  with  6  or  15 
MV  photon  energy.  A  wedge  filter  was  used  in  the  lateral  direction.  When  a  physical  wedge  was 
used,  we  avoided  placing  it  in  the  medial  field  to  reduce  the  scatter  dose  to  the  lung/heart  and  the 
contra-lateral  breast.  In  this  case,  both  fields  were  modulated  with  multiple  segments.  The 
segmented  fields  in  the  lateral  direction  were  delivered  concurrently  with  the  physical  wedge  in 
place.  All  plans  were  obtained  through  manual  trial-and-error  process.  The  MART  planning  was 
done  with  Corvus  system  (North  American  Scientific,  Cranberry  Township,  PA).  Several  closely 
related  clinical  issues,  including  the  type  of  wedges  (physical  or  dynamic),  field  size, 
incorporation  of  MLC  transmission  into  the  step-and-shoot  delivery,  and  QA,  were  extensively 
investigated  during  the  implementation  of  the  MART  technique.  Our  study  indicated  that  the 
MART  markedly  improves  breast  irradiation  and  provides  superior  dose  distributions  needed  to 
reduce  radiation  side  effects  and  complications.  The  technique  is  especially  valuable  for 
radiation  treatment  of  large-breasted  women,  where  it  is  difficult  to  achieve  homogeneous  target 
dose  distribution. 

Fig.  2  showed  a  typical 
clinical  case  that  fell  into  the 
category  of  intermediately 
complicated  or  complicated  cases. 

In  this  comparison  study,  the 
prescribed  dose  was  specified  to  a 
point  ~3  cm  anterior  to  the 
isocenter  and  it  was  desired  that  the 
100%  isodose  curve  to  cover  the 
breast  target  volume.  The 
prescription  dose  of  all  treatment 
plans  were  scaled  to  5,040  cGy. 

The  hot  spots  of  the  standard  OTF 
plans  in  the  breast  volume  ranged 
from  109%  to  118%  when 
normalized  to  the  prescribed  dose. 

As  shown  in  Fig.  1,  the  isodose 

distributions  in  the  central 

transverse  section  and  in  a  plane 

perpendicular  to  the  incident  beams 
for  the  standard  OTF  and  MART  treatments,  indicate  that  the  dose  inhomogeneity  in  the  target 

volume  was  significantly  reduced  with  MART,  as  well  as  reduction  in  the  high  dose  to  the 

ipsilateral  lung  and  heart  when  compared  with  the  OTF  plan.  The  target  maximum  dose  was 

reduced  from  118%  to  1 12%  for  the  MART  plan.  Furthermore,  the  target  volume  receiving  high 

dose  irradiation  was  significantly  reduced.  In  order  to  include  the  medial  breast  tissue  into  the 

radiation  field,  ~10%  of  the  heart  volume  and  the  left  lung  were  included  in  the  tangential  fields. 

As  thus,  a  significant  fraction  of  the  heart  and  lung  receives  high  radiation  dose  in  standard  OTF 

plan.  The  high  doses  to  the  heart  was  reduced  by  almost  6%  using  MART  technique  as  a  result 

of  adding  one  additional  segment  in  each  incident  beam,  together  with  ~5%  improvement  in  the 

maximum  target  dose  in  the  target  volume.  The  heart  volume  and  ipsilateral  lung  volumes 


(a)  (b) 


Figure  2  Comparison  of  the  isodose  distributions  of  the  treatment 
plans  of  the  left-sided  breast  case  using  the  tangential  field 
technique  (a),  MART  (b).  Target  volume  includes  the  whole  breast 
and  the  internal  mammary  nodes.  Isodose  levels  are  shown  at 
1 10%,  100%,  90%,  70%,  50%,  30%,  and  10%. 
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receiving  high  dose  irradiation  were  also  markedly  reduced  for  MART  treatment.  In  addition  to 
the  comparison  of  MART  with  OTF,  the  IMRT  plans  were  subsequently  made  for  all  the  35 
patients.  Generally,  the  dose  distributions  of  IMRT  is  similar  to  MART  plans,  however,  the  time 
cost  of  MART  is  significantly  less. 


D.  Efficiency  and  accuracy  of  MART  ( Nov  2005  ~  Nov  2006) 

For  the  MART  assessment,  we  have  conducted  many  numerical  and  physical  phantom 
experiments,  especially  focusing  on  the  breathing  motion  problems.  In  radiation  therapy, 
respiratory  motion  poses  significant  challenges  for  tumors  in  breast.  The  motion  can  distort  the 
shape  of  an  object,  degrade  the  anatomic  position  reproducibility  during  imaging,  and  necessitate 
larger  margins  during  radiotherapy  planning.  It  also  causes  inaccuracy  in  estimating  the  tumor 
volume,  thereby  preventing  an  effective  dose  escalation  for  the  treatment  of  a  target  tumor. 
These  issues  make  it  difficult  to  achieve  the  desired  goals  of  conformal  radiotherapy. 

We  have  implemented  a  gated 
irradiation  technique  for  the  treatment 
of  left-sided  breast  cancer  at  Stanford 
Hospital,  where  the  goal  is  to  better 
conform  the  tumor  target  while  sparing 
the  heart.  In  figure  3  we  show  an 
example  of  the  treatment.  Fig.  3a  and 
3b  show  the  inhale  and  exhale  phases, 
respectively.  Fig.  3  c  shows  the  plan  for 
gated  treatment  of  the  patient.  As  seen 
from  the  isodose  plot,  the  heart  of  the 
patient  is  spared  greatly  as  compared 
with  conventional  free-breathing 
treatment.  Figure  3d  shows  the 
patient’s  breathing  pattern.  A  temporal 
margin  of  20%  (the  portion  highlighted 
by  thicker  lines  on  the  breathing  curve) 
was  imposed  during  the  gated  radiation 
delivery  process. 

The  above  gated  treatment 
requires  a  Four-dimensional  (4D) 
simulation  CT.  4D  CT  scans,  acquired 
synchronously  with  a  respiratory  signal,  provide  not  only  the  three-dimensional  (3D)  spatial 
information,  but  also  temporal  changes  of  the  anatomy  as  a  function  of  the  respiratory  phase 
during  the  imaging,  and  can  therefore  be  employed  in  4D  treatment  planning  to  explicitly 
account  for  the  respiratory  motion,  for  example,  the  respiration  gated  IMRT/MART  treatment, 
where  the  photon  beams  are  only  switched  on  in  particular  phase  window  determined  by  a  real¬ 
time  motion  tracking  system  (RPM).  An  astonishing  problem  with  4D  CT  for  breast  cancer 
patients  is  the  extremely  high  radiation  dose  during  the  4D  imaging  (about  10  to  20  times  higher 
than  a  regular  3D  CT  scan).  In  particular  for  patients  with  partial  breast  treatment,  such  high  dose 
screening  may  lead  to  secondary  radiation-induced  cancer  in  the  normal  tissue  (the  other  normal 
breast  for  example).  This  concern  led  to  our  extra  contribution  beyond  the  initial  aims  of  this 
project:  a  novel  technique  to  lower  the  radiation  exposure  to  the  patient  in  4D  imaging  while 


Figure  3  Gated  breast  radiation  treatment  with  4D  CT.  (a)  and  (b) 
are  inspiration  and  expiration  phases,  respectively,  (c)  is  the  gated 
treatment  plan,  and  (d)  shows  the  patient  breathing  pattern,  where 
the  thick  green  line  indicates  the  20%  to  80%  respiratory  phase 
window. 
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maintaining  the  advantage  of  the  provided  kinetic  information  of  the  tumor  and  organs  at  risk 
(OAR)  by  4D  CT.  The  method  that  we  have  developed  for  this  purpose  is  to  perform  4D  CT 
scans  at  relatively  low  x-ray  tube  current,  hence  reducing  the  radiation  exposure  of  the  patients. 
To  deal  with  the  increased  statistical  noise  caused  by  the  low  current,  a  novel  4D  penalized 
weighted  least  square  (4D-PWLS)  smoothing  method  has  been  proposed,  which  can  incorporate 
both  spatial  and  phase  information.  The  4D  images  at  different  phases  were  registered  to  the 
same  phase  via  a  deformable  model,  and  a  regularization  term  combining  temporal  and  spatial 
neighbors  can  be  designed  for  the  4D-PWLS  objective  function. 

The  proposed  method  was  validated  with  phantom  experiments  and  applied  to  patient 
studies,  where  superior  noise  suppression  and  resolution  preservation  were  observed.  4D  CT 
images  were  acquired  with  a  combined  PET/CT  scanner  (Discovery  ST/LightSpeed  8-slice, 
General  Electric  Medical  Systems)  in  our  clinic.  Fig.  4  show  the  phantoms  used  for  the 
validation  of  our  4D-PWLS  method:  one  is  a  commercial  calibration  phantom  CatPhan®  600 
(The  Phantom  Laboratory,  Inc.,  Salem,  NY),  and  the  other  is  an  anthropomorphic  thorax 
phantom.  To  acquire  the  4D  CT  data,  each  phantom  was  placed  on  the  top  of  a  platform  capable 
of  sinusoidal  motion  along  the  cranial-caudal  direction.  A  Real-time  Position  Management 
(RPM)  respiratory  gating  system  (Varian  Medical  Systems,  Palo  Alto,  CA)  was  used  to  record 
the  motion  by  tracking  two  infrared  reflective  markers,  rigidly  mounted  on  a  plastic  block  on  the 
top  of  the  phantom,  by  means  of  an  infrared  video  camera  mounted  on  the  PET/CT  table.  The 
clinical  4D-CT  patient  studies  were  performed  on  the  same  scanner.  Patients  were  asked  to 
breathe  normally  during  the  scan.  The  plastic  block  with  two  infrared  reflective  markers  was 
taped  on  the  top  of  the  patient's  abdomen,  placed  medially  and  a  few  cm  inferior  to  the  xiphoid 
processes.  The  respiratory  signal  of  the  patient  from  the  RPM  was  recorded  and  synchronized 
with  CT  data  acquisition,  similar  to  the  phantom  scan.  After  the  scan  data  were  prospectively 
reconstructed  at  the  PET/CT  scanner,  both  the  CT  images  and  the  corresponding  motion  data 
recorded  by  the  RPM  system  were  transferred  to  a  GE  Advantage  Workstation  (GE  Medical 
Systems,  Waukesha,  WI).  The  "Advantage  4D"  software  on  the  workstation  simultaneously 
displays  the  CT  images  and  the  motion  data,  and  sorts  the  cine  images  into  a  set  of  respiratory 
phase  images.  The  CatPhan  study  showed  that  the  new  method  improved  the  signal-to-noise 
(SNR)  about  2.5x,  with  only  a  6%  loss  in  spatial  resolution.  In  the  human  thorax  phantom  study, 
The  effective  reduction  of  the  image  noise  are  also  observed,  where  the  average  SNR  of  10  mA 
images  increased  by  more  than  three-fold  from  0.051  to  0.165  after  the  proposed  post-processing. 
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Figure  4.  Phantoms  used  for  the  motion  study  (a)  CatPhan,  (b)  anthropomorphic  thorax  phantom  and  the  motion 
platform. 
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One  example  of  patient  studies  is  shown  in  Fig  5.  The  left  column  is  the  original  phase 
image  (end-inspiration  phase)  obtained  from  GE  Advantage  Workstation,  and  the  right  column  is 
the  image  after  processing  with  the  proposed  4D-PWLS  enhancement.  Successful  noise 
suppression  is  observed.  For  this  patient,  the  SNRs  increased  from  2.204  to  4.558  for  the  end- 
expiration  phase,  and  from  1.741  to  3.862  for  the  end-inspiration  phase  in  the  selected  ROIs. 


Figure  5.  One  example  of  the  patient  studies  for  the  4D-PWLS  method.  The  phase  shown  here  is  the  end-inspiration 
phase.  The  left  is  the  regular  image  obtained  from  the  GE  4D  Advantage  Workstation,  and  the  right  shows  the  image 
after  4D-PWLS  processing.  The  red  rectangles  represent  the  selected  ROI  for  the  calculation  of  SNRs,  each  of  which 
contain  5X5X5  voxels,  where  the  SNR  is  improved  from  2.2  to  4.6. 


KEY  RESEARCH  ACCOMPLISHMENTS 


•  Developed  and  clinically  implemented  a  manual  MART  planning  procedure  at  Stanford 
Hospital. 

•  Developed  an  aperture -based  optimization  algorithm  for  IMRT/MART  treatment  planning, 
and  assessed  the  dosimetric  improvement  resulted  from  the  optimized  MART. 

•  Improved  the  general  inverse  planning  framework  by  introducing  a  voxel  dependent 
penalty  and  developed  a  dedicated  genetic  algorithm  for  MART  planning.  Demonstrated  its 
superiority  over  conventional  treatment  planning  systems  in  terms  of  dose  homogeneity 
and  efficiency. 

•  Compared  35  breast  cancer  patient  cases  with  both  MART  and  traditional  OTF  plans,  and 
demonstrated  the  advantage  of  MART  over  OTF.  Now  over  100  breast  cancer  patients 
have  been  treated  using  the  new  technique. 

•  Developed  a  novel  strategy  of  gated  IMRT/MART  treatment  for  breast  caner.  Constructed 
a  breathing  motion  phantom  for  the  4D  imaging/treatment  research. 

•  Developed  a  method  that  significantly  reduced  the  radiation  exposure  to  the  patient  during 
the  4D  CT  imaging,  which  is  an  essential  step  in  4D  radiation  therapy  for  breast  cancer. 

•  Developed  a  method  that  significantly  improved  the  4D  PET  imaging  quality. 
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•  Li  T,  Xing  L,  Munro  P,  McGuinness  C,  Chao  M,  Yang  Y,  Loo  B,  and  Koong  A,  “Four¬ 
dimensional  cone-beam  computed  tomography  using  an  on-board  imager,”  Medical  Physics 
33(10),  3825-3833,  2006. 

•  Li  T,  Thomdyke  B,  Schreibmann  E,  Yang  Y,  and  Xing  L,  “Model-based  image  reconstmction  for 
four-dimensional  PET,”  Medical  Physics  33(5),  1288-1298,  2006. 

•  Li  T,  Schreibmann,  E,  Yang  Y  and  Xing  L,  “Motion  correction  for  improved  target  localization 
with  on-board  cone-beam  computed  tomography,”  Physics  in  Medicine  Biology  51(2),  253-267, 
2006. 

•  Li  T,  Schreibmann  E,  Thomdyke  B,  Tillman  G,  Boyer  A,  Koong  A,  Goodman  K  and  Xing  L, 
“Radiation  dose  reduction  in  4D  computed  tomography,”  Medical  Physics  32(12),  3650-3660, 
2005. 

•  Shou  Z,  Yang  Y,  Cotmtz  C  and  Xing  L,  “Quantitation  of  the  a  priori  dosimetric  capabilities  of 
spatial  points  in  inverse  planning  and  its  significant  implication  in  defining  IMRT  solution 
space,”  Physics  in  Medicine  and  Biology  50(7),  1469-1482,  2005. 

Published  Abstracts: 

The  Pis’  group  has  also  been  active  in  disseminating  our  research  results.  The  following  are 

some  of  the  presentations  given  in  various  national/intemational  meetings. 

•  Li  T,  Cotmtz  C,  Gofinett  D,  and  Xing  L,  “Segmenation-based  breast  IMRT  using  a  genetic  dose 
optimization  algorithm,”  Era  of  Hope  (Department  of  Defense  Breast  Cancer  Research  Program 
Meeting)  131,  2005. 

•  Yang  Y,  Li  T,  Schreibmann,  Boyer  A,  and  Xing  L,  “Is  cone  beam  CT  suitable  for  dose 
verification?”  Medical  Physics  32,  2160,  2005. 

•  Li  T,  Yang  Y,  Schreibmann,  and  Xing  L,  “A  new  cone-beam  CT  repositioning  technique  through 
deformable  registration,”  Medical  Physics  32,  2160,  2005. 

•  Xing  L,  Schreibmann,  Yang  Y,  Boyer  A,  and  Li  T,  “Image  segmenatation  in  4D  CT  based  on  a 
deformable  image  registration  model,”  Medical  Physics  32,  2095,  2005. 

•  Kim  G,  Li  T,  Yang  Y,  Schreibmann,  Thomdyke  B,  Boyer  A,  and  Xing  L,  “Influence  of 
respiratory  motion  on  cone-beam  CT  imaging  of  thorax  and  abdomen,”  Medical  Physics  32, 
1934,2005. 

•  Thomdyke  B,  Schreibmann,  Li  T,  Boyer  A,  and  Xing  L,  “A  comparison  of  amplitude-  and  phase- 
based  4D  CT,”  Medical  Physics  32,  1919,  2005. 

•  Lo  A,  Shou  Z  and  Xing  L,  Quantitative  comparison  of  aperture-based  and  beamlet-based  inverse 
planning  techniques,  International  Journal  of  Radiation  Oncology,  Biology,  Physics,  Volume  60, 
1,  Supplement  1,  September  2004,  Pages  S630-S630; 

•  Shou  Z  and  Xing  L,  Aperture-based  IMRT  inverse  planning  with  incorporation  of  organ  motion, 
International  Journal  of  Radiation  Oncology,  Biology,  Physics,  Volume  60,  1,  Supplement  1, 
September  2004,  Pages  S631-S632. 

•  Xing  L  and  Shou  Z,  Intrinsic  spatial  heterogeneity  in  inverse  planning  and  its  significant  role  in 
defining  the  universe  of  IMRT  solution  space,  International  Journal  of  Radiation  Oncology, 
Biology,  Physics,  Volume  60,  1,  Supplement  1,  September  2004,  Pages  S635-S635. 

•  Shou  Z  and  Xing  L:  Improve  IMRT  distribution  by  using  spatially  non-uniform  important 
factors,  oral  presentation  in  2004  annual  meeting  of  ICCR  in  Seoul,  Korea.  120-123,  2004. 
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•  Xing  L,  Hunjan  S,  Lian  J,  Yang  Y,  Shou  Z,  Schreibmann  E,  Boyer  A,  Towards  Biologically 
Conformal  Radiation  Therapy:  Functional  and  Molecular  Image  Guided  Intensity  Modulated 
Radiation  Therapy,  XIVth  International  Conference  on  the  Use  of  Computers  in  Radiation 
Therapy  (ICCR),  Soul,  Korea,  2004. 

CONCLUSIONS 

Compared  with  conventional  methods  for  breast  cancer  RT  such  as  OTF,  modem  IMRT 
provides  unprecedented  potential  to  deliver  3D-dose  distributions  with  superior  tumor 
conformality  and  normal  tissue  spare.  In  practice,  however,  the  current  IMRT  process  deviates 
significantly  from  the  conventional  approach  and  requires  considerable  costs  of  the  time  and 
health  care.  The  complexity  of  the  treatment  and  deliver  procedure  make  this  technology 
problematic  for  the  breast  radiotherapy.  We  have  carried  out  a  systematic  study  on  a  new 
solution,  which  is  a  hybrid  of  the  advantages  of  the  IMRT  and  OTF,  termed  as  multiple  aperture 
radiotherapy  (MART).  The  efficiency  in  dose  shaping  and  simplicity  in  dose  delivery  make  it  a 
perfect  choice  for  breast  cancer  treatment.  A  number  of  important  milestones  have  been 
accomplished,  which  include  (i)  developed  an  efficient  optimization  method  for  MART  inverse 
planning;  (ii)  implemented  MART  in  a  clinic  environment  for  breast  treatment;  (iii)  proposed  a 
gated  IMRT/MART  based  on  the  cutting-edge  technology  of  4D  CT  imaging;  (iv)  established  a 
robust  technique  for  low-dose  4D  CT  acquisition  to  reduce  the  risks  of  radiation-induced  cancer. 

It  is  expected  these  tools  will  greatly  facilitate  the  imaging,  planning,  delivery,  and  quality 
assurance  of  breast  radiation  treatment. 
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On-board  cone-beam  computed  tomography  (CBCT)  has  recently  become  available  to  provide 
volumetric  information  of  a  patient  in  the  treatment  position,  and  holds  promises  for  improved 
target  localization  and  irradiation  dose  verification.  The  design  of  currently  available  on-board 
CBCT,  however,  is  far  from  optimal.  Its  quality  is  adversely  influenced  by  many  factors,  such  as 
scatter,  beam  hardening,  and  intra- scanning  organ  motion.  In  this  work  we  quantitatively  study  the 
influence  of  organ  motion  on  CBCT  imaging  and  investigate  a  strategy  to  acquire  high  quality 
phase-resolved  [four-dimensional  (4D)]  CBCT  images  based  on  phase  binning  of  the  CBCT  pro¬ 
jection  data.  An  efficient  and  robust  method  for  binning  CBCT  data  according  to  the  patient’s 
respiratory  phase  derived  in  the  projection  space  was  developed.  The  phase-binned  projections  were 
reconstructed  using  the  conventional  Feldkamp  algorithm  to  yield  4D  CBCT  images.  Both  phantom 
and  patient  studies  were  carried  out  to  validate  the  technique  and  to  optimize  the  4D  CBCT  data 
acquisition  protocol.  Several  factors  that  are  important  to  the  clinical  implementation  of  the  tech¬ 
nique,  such  as  the  image  quality,  scanning  time,  number  of  projections,  and  radiation  dose,  were 
analyzed  for  various  scanning  schemes.  The  general  references  drawn  from  this  study  are:  (i) 
reliable  phase  binning  of  CBCT  projections  is  accomplishable  with  the  aid  of  external  or  internal 
marker  and  simple  analysis  of  its  trace  in  the  projection  space,  and  (ii)  artifact-free  4D  CBCT 
images  can  be  obtained  without  increasing  the  patient  radiation  dose  as  compared  to  the  current  3D 
CBCT  scan.  ©  2006  American  Association  of  Physicists  in  Medicine.  [DOI:  10.1118/1.2349692] 

Key  words:  cone-beam,  4D  CT,  on-board  imager,  IGRT,  organ  motion 


I.  INTRODUCTION 

Modern  radiation  therapy  techniques,  such  as  three- 
dimensional  (3D)  conformal  radiotherapy  and  intensity- 
modulated  radiation  therapy  (IMRT),  provide  unprecedented 
means  for  producing  exquisitely  shaped  radiation  doses  that 
closely  conform  to  the  tumor  dimensions  while  sparing  sen¬ 
sitive  structures.  As  a  result  of  greatly  enhanced  dose  confor¬ 
mality,  more  accurate  beam  targeting  becomes  an  urgent  is¬ 
sue  in  radiation  therapy.  The  need  to  improve  targeting  in 
high  precision  radiation  therapy  has  recently  spurred  a  flood 
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of  research  activities  in  image-guided  radiation  therapy. 
Cone-beam  computed  tomography  (CBCT),  based  upon  flat- 
panel  technology,  integrated  with  a  medical  linear  accelerator 
has  recently  become  available  from  Linac  vendors  for 
therapy  guidance.  The  volumetric  images  provided  by  CBCT 
are  being  used  to  verify  and  correct  the  planning  patient 
setup  in  the  Linac  coordinates  by  comparing  with  the  patient 
position  defined  in  the  treatment  plan.8-10  In  addition,  CBCT 
data  acquired  prior  to  the  treatment  can  be  used  to  recalculate 
and  verify  the  delivered  dose  to  the  patient  on  that  treatment 
day,  or  might  even  be  used  to  generate  online  treatment 


plans.11-14  All  these  applications  strongly  depend  on  the 
quality  of  the  CBCT  images,  which  is  often  severely  affected 
by  artifacts  caused  by  intra- scanning  organ  motion.  Accord¬ 
ing  to  the  current  IEC  standard,  on-board  CBCT  has  a  lim¬ 
ited  rotation  speed  about  60  s  per  rotation.  A  full  (360°)  scan 
therefore  consists  of  projection  data  from  about  20  respira¬ 
tory  cycles  of  the  patient,  resulting  in  large  amount  of  incon¬ 
sistency  in  the  CBCT  projection  data.  Reconstruction  based 
on  theories  that  assume  stationary  object  thus  generates  lots 
of  artifacts,  which  can  be  more  severe  than  the  conventional 
diagnostic  CT.  The  artifacts  not  only  make  it  difficult  to  see 
the  extent  of  the  tumor,  but  also  inhibit  direct  use  of  CBCT 
for  dose  calculation.14,15 

Methods  to  account  for  respiratory  motion  during  CT  im¬ 
aging  can  be  divided  into  two  categories:  one  is  to  find  new 
reconstruction  theory  that  can  model  and  compensate  for  the 
motion  effects,16-19  and  the  other  is  to  seek  new  acquisition 
procedures  to  relieve  the  motion  effects.  The  later  is  repre¬ 
sented  by  techniques  such  as  breath-hold,  respiratory  gating 
and  four-dimensional  (4D)  CT.  Implementation  of  4D  CT  for 
conventional  CT  scanner  has  been  successfully  demonstrated 
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recently,  for  example  with  cine-mode  protocol, 20-23  but  there 
has  been  limited  literatures  on  discussing  4D  acquisition 
with  on-board  CBCT  imager.24,25  The  major  difference  here 
is  that,  with  conventional  diagnostic  CT,  each  single-rotation 
scan  can  be  completed  in  a  time  interval  (^0.5  s),  much 
shorter  than  a  typical  breathing  period.  Therefore,  the  projec¬ 
tions  obtained  in  each  individual  rotation  have  “complete 
data”  for  an  artifact-free  reconstruction,  and  reconstructed 
images  from  multiple  scans  can  be  subsequently  sorted  to 
generate  a  4D  sequence.  For  on-board  CBCT,  however,  pro¬ 
jection  data  of  each  single  rotation  have  been  contaminated 
with  motion  already,  and  no  “pure-phase”  image  can  be  de¬ 
rived  directly.  A  recently  proposed  method  of  the  “respiration 
correlated  CBCT”  sorts  the  projection  data  before  image 
reconstruction24  and  provides  trajectories  of  the  involved 
structures.  However,  strong  view  aliasing  artifacts  are 
present  in  so  obtained  4D  CBCT  images  due  to  the  reduced 
number  of  projections  in  each  phase. 

The  aim  of  this  research  is  to  investigate  an  effective 
strategy  of  acquiring  high  quality  4D  CBCT  scans  to 
eliminate/reduce  the  respiration-induced  artifacts  while  ob¬ 
taining  temporal  information  of  the  patient  anatomy.  In  the 
following,  we  describe  detailed  step-by-step  procedures  of 
data  acquisition,  respiration  signal  extraction,  and  image  re¬ 
construction.  Different  scan  settings  have  been  studied  for  a 
motion  phantom  by  varying  the  x-ray  tube  currents  and  the 
number  of  revolutions.  Quantitative  comparisons  of  the  data 
are  performed,  and  the  clinically  most  feasible  4D  scan  set¬ 
tings  are  explored.  Finally,  4D  CBCT  study  of  a  lung  cancer 
patient  is  presented. 

II.  METHODS  AND  MATERIALS 
A.  4D  CBCT  data  acquisition 

The  on-board  CBCT  imaging  system  used  in  this  work 
was  the  ExactArms  (kV  source/detector  arms)  of  a  Trilogy™ 
treatment  system  (Varian  Medical  Systems,  Palo  Alto,  CA). 
The  current  clinical  protocol  for  acquiring  a  regular  3D 
CBCT  scan  is  with  125  kVp  voltage  and  80  mA  current.  The 
gantry  rotation  speed  was  6  °  / s,  and  the  duration  of  the  x-ray 
pulse  at  each  projection  angle  was  25  ms.  A  full  rotation 
(slightly  over  360°)  consisted  of  over  680  projections  with  an 
angle  interval  of  about  0.53°.  The  dimension  of  each  ac¬ 
quired  projection  image  was  about  397  mm  X  298  mm,  con¬ 
taining  1024  X  768  pixels.  The  system  has  a  field  of  view 
(FOV)  of  25  cm  X  25  cm  in  the  transverse  plane  and  17  cm 
in  the  longitudinal  direction  (full-fan  mode),  which  can  be 
increased  to  50  cm  X  50  cm  in  the  transverse  plane  by  shift¬ 
ing  the  detector  laterally  (half-fan  mode).  The  acquisition 
time  of  one  rotation  was  about  1.07  min,  and  the  radiation 
dose  of  such  acquired  CBCT  was  around  3.8  cGy  for  body 
scan  and  9.0  cGy  for  head  scan. 

Since  a  full  CBCT  scan  usually  spans  several  respiration 
cycles,  one  can  sort  the  collected  projections  into  a  few 
groups  according  to  their  respiration  phases  and  then  recon¬ 
struct  each  group  of  the  “phase-binned”  projection  data  to 
derive  4D  CBCT  images.  One  problem  in  doing  so  is  that 
each  phase  group  contains  much  fewer  projections  than  a 


regular  scan,  which  causes  reconstruction  artifacts  due  to  in¬ 
sufficient  angular  sampling  and  dramatically  degrades  the 
image  quality.  Crucial  issues  in  developing  4D  CBCT  are 
indeed  how  to  enhance  the  amount  of  information  in  each 
phase  and  how  to  deal  with  the  limited  number  of  projections 
in  each  phase.  Here  we  investigate  an  intuitive  technique  of 
increasing  the  projection  number  by  using  multiple-rotation 
scan.  But  the  methodology  developed  in  the  following  is 
well  suited  for  related  approaches  based  on,  for  instance, 
“slowing  down  the  gantry  rotation,”  which  makes  a  full  scan 
contain  data  over  more  respiration  cycles. 

B.  Extract  respiration  signal,  phase  sorting,  and 
reconstruction 

To  sort  the  acquired  CBCT  data  in  projection  space,  each 
projection  needs  to  be  stamped  with  a  respiration  phase  dur¬ 
ing  the  data  acquisition.  This  can  be  done  by  synchronizing  a 
real-time  motion  tracking  system  (for  example  the  RPM) 
with  the  CBCT  acquisition,  similar  to  that  in  the  conven¬ 
tional  4D  CT  acquisitions.21,22  In  this  work,  we  used  another 
simple  and  reliable  method  to  achieve  this  by  monitoring  the 
trajectory  of  a  radiopaque  fiducial  marker  (radius  of  1  mm) 
adhered  to  the  patient  skin  in  the  projection  space.  It  is  nec¬ 
essary  to  place  the  fiducial  marker  to  where  it  can  be  seen  by 
the  detector  at  any  angle.  In  the  case  that  a  patient  size  is  too 
large  and  the  marker  exceeds  the  FOV  at  certain  angles,  a 
second  marker  placed  at  a  different  position  can  be  used. 
Many  studies  have  shown  that  the  skin  marker  describes  the 
respiration  movement  well  and  can  serve  as  a  surrogate  char¬ 
acterizing  the  patient’s  breathing  motion.26,27 

To  automatically  identify  the  marker  in  the  projection  im¬ 
ages,  we  have  developed  a  fast  and  robust  “successive 
searching”  process.  The  analysis  starts  from  the  first  projec¬ 
tion  with  an  initial  guess  of  the  marker  position,  which  is 
calculated  based  on  the  cone-beam  imaging  system  geometry 
(the  kV  x-ray  source/flat-panel  detector  positions),  and  the 
approximate  marker  coordinates  in  the  treatment  reference 
frame.  Furthermore,  a  region-of-interest  (ROI)  window  of 
size  of  40  X  40  pixels  is  applied  around  the  estimated  marker 
position  in  the  projection  image  to  ensure  the  maker  is  in¬ 
cluded.  The  exact  marker  location  is  then  determined  by 
searching  for  the  minimal  intensity  in  the  window,  and  is  set 
to  be  the  center  of  the  ROI  window  for  the  next  projection 
because  of  the  proximity  of  the  two  successive  projections. 
Once  the  ROI  is  placed,  the  searching  process  is  repeated  to 
find  the  exact  marker  location  in  the  new  projection.  Al¬ 
though  the  successively  updated  window  generally  follows 
the  marker  from  one  projection  to  another,  due  to  the  pres¬ 
ence  of  dense  objects,  it  is  possible  that  the  minimal  intensity 
in  the  selected  window  does  not  correspond  to  the  marker 
location  in  certain  situation.  To  make  the  searching  process 
more  robust,  we  used  a  blurred  mask  subtraction  for  an  edge 
enhancement.  Typical  time  for  finding  the  marker  in  a  pro¬ 
jection  image  is  less  than  0.1  s  using  matlab  (Pentium  4 
CPU  2.66  GHz,  512  MB  RAM  PC  platform.)  The  error  of 
the  marker  location  would  be  1-2  pixels  with  the  detector 
pixel  size  of  0.388  mm.  It  should  be  noted  that  the  error  was 
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Fig.  1 .  CatPhan®  600  phantom  and  its  module  specifi¬ 
cation.  The  module  section  CTP404  was  used  in  this 
work  to  demonstrate  the  image  quality  visually. 


measured  only  for  the  periodically  moving  phantom,  which 
represented  a  best-case  scenario. 

After  the  marker  is  identified  in  the  projection  images,  the 
phase  of  each  projection  is  then  determined  according  to  the 
marker  location.  In  this  work,  the  end  inspiration  is  defined 
as  the  0%  phase  with  100%  corresponding  to  a  complete 
breathing  cycle.  The  aim  of  phase  binning  is  to  group  pro¬ 
jections  of  similar  phase  together.  For  any  specific  phase,  a 
phase  uncertainty  8  is  user  defined  to  control  the  phase  range 
or  residual  motion  contained  in  the  phase-binned  data.  After 
phase  binning,  the  data  can  be  reconstructed  with  the  con¬ 
ventional  Feldkamp  algorithm  to  get  the  4D  images.28  In  this 
work,  the  size  of  each  reconstructed  volumetric  image  was 
512X512X161,  with  voxel  size  of  0.5  X  0.5  X  1.0  mm3.  It 
should  be  noted  that,  the  sorted  projections  usually  distribute 
unevenly  around  the  360°  circle,  a  correct  weight  should  be 
used  for  each  projection  according  to  the  gantry  angle  incre¬ 
ment  between  projections  to  account  for  the  nonuniformity. 


C.  Phantom  study 

A  series  of  experiments  using  a  motion  phantom  were 
performed  to  validate  the  proposed  4D  CBCT  technique  and 
to  search  for  the  optimal  data  acquisition  protocol.  The  mo¬ 
tion  phantom  consists  of  a  commercial  CT  calibration  phan¬ 
tom  CatPhan®  600  (The  Phantom  Laboratory,  Inc.,  Salem, 
NY)  placed  on  the  top  of  a  platform  capable  of  sinusoidal 
motion  along  three  directions  (see  Fig.  1).  A  small  metal 
marker  of  1  mm  size  was  attached  to  the  top  of  the  phantom 
during  the  scans  as  an  indicator  of  the  motion  status  in  the 
projection  images.  Two  motion  modes  of  the  phantom  sys¬ 
tem  were  employed  to  evaluate  the  proposed  4D  CBCT  ac¬ 
quisition  technique.  In  the  first  one,  the  phantom  was  set  in 
motion  in  the  longitudinal  direction  [superior-inferior  (SI) 
direction]  with  an  amplitude  of  1.0  cm  (or  peak-valley  dis¬ 
placement  of  2.0  cm).  The  motion  period  was  set  to  5.5  s.  In 
the  second  mode,  the  phantom  moved  with  the  maximal  dis¬ 


placements  5.5  cm  in  the  SI  direction,  1.5  cm  in  the  anterior- 
posterior  direction,  and  0.5  cm  in  the  lateral  direction.  The 
moving  period  was  4.6  s. 

The  first  motion  mode  was  used  to  study  the  dependence 
of  4D  CBCT  quality  on  the  number  of  gantry  rotations.  For 
this  purpose,  up  to  eight  CBCT  scans  of  the  motion  phantom 
were  acquired  with  the  gantry  rotated  clockwise  and  coun¬ 
terclockwise  alternatively.  The  voltage  and  x-ray  pulse  width 
were  kept  the  same  at  120  kVp  and  25  ms,  respectively.  The 
x-ray  tube  current  was  set  to  be  10  mA  for  each  of  the  eight 
scans.  The  4D  CBCT  images  using  N  rotations  were  then 
reconstructed  based  on  the  projection  sorting  scheme  de¬ 
scribed  earlier,  where  N=  1,2, 3,...  ,8.  For  comparison,  3D 
stationary  CBCT  scans  were  also  performed  at  two  phantom 
positions  corresponding  to  the  “mid-expiration”  and  “end- 
expiration”  phases.  The  quality  of  the  4D  CBCT  images  cor¬ 
responding  to  different  number  of  rotations  was  analyzed 
quantitatively  using  a  relative  error  metric,  which  will  be 
described  later. 

In  the  second  set  of  experiments,  a  number  of  4D  acqui¬ 
sitions  were  carried  out  with  the  radiation  dose  kept  at  the 
same  level.  Specifically,  scans  of  one  rotation  with  80  mA, 
two  rotation  with  40  mA,  four  rotation  with  20  mA,  and 
eight  rotation  with  10  mA  were  acquired.  Again,  for  quanti¬ 
tative  comparison,  two  additional  scans,  with  the  phantom 
“frozen”  at  corresponding  positions  of  the  mid-expiration 
and  end-expiration  phases,  were  obtained  using  80  mA  cur¬ 
rent. 

D.  Patient  data  acquisition 

Four-dimensional  CBCT  scan  was  performed  for  a  lung 
cancer  patient  who  underwent  gated  IMRT  treatment  in  our 
clinic.  In  this  acquisition,  to  ensure  the  quality  of  the  imag¬ 
ing  study,  four  continuous  gantry  rotations  were  done  with 
the  tube  current  set  to  32  mA.  The  full-fan  scan  mode  of  the 
Trilogy  was  employed  for  the  study,  in  which  the  maximum 
FOV  in  transverse  plane  was  the  25  cmX  25  cm,  and  17  cm 
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in  longitudinal  direction.  A  full-fan  (symmetric)  bow-tie  fil¬ 
ter  was  used  between  the  kV  x-ray  source  and  the  patient. 
The  average  breathing  cycle  of  the  patient  was  4.2  s.  The  4D 
CBCT  scan  added  about  5  min  to  the  treatment  time.  After 
the  data  acquisition,  the  projections  were  sorted  according  to 
their  phases  and  subsequently  reconstructed  into  four  single 
phases,  namely,  the  end-inspiration  (phase  0%),  mid¬ 
expiration  (phase  25%),  end-expiration  (phase  50%),  and 
mid-inspiration  phases  (phase  75%). 

E.  Image  quality  assessment  metric 

Since  the  images  in  this  work  involve  significant  artifacts, 
common  metric  such  as  contrast- to-noise  ratio  (CNR),  de¬ 
fined  as  CNR=|(51-52)|/o',  may  not  be  suitable,  in  the  sense 
that  a  high  CNR  output  may  not  indicate  a  better  image.  For 
example,  dark  streak  artifacts  may  accidentally  decrease  the 
mean  value  of  S2,  leading  to  a  high  CNR  value.  To  quantita¬ 
tively  assess  the  images  obtained  from  different  CBCT  set¬ 
tings,  we  chose  the  “relative  root-mean-square  error”  (de¬ 
noted  by  RE)  as  a  figure  of  merit  of  the  image  quality,  which 
is  defined  as 

RE  =  ^2  [WO  -  S80  mA«]2  j  2  [Sm(0  -  S80  mA(0]2 j  1/2, 

where  S4D  denotes  the  4D  single-phase  image,  SM  and  S8 0  mA 
are  the  standard  80  mA  3D  CBCT  images  of  the  phantom 
with  and  without  motion,  respectively,  and  the  summation 
runs  over  all  voxels.  In  the  above  equation,  580  mA  image  is 
used  as  the  “gold  standard,”  and  the  mean  square  error  of 
between  the  4D  images  and  the  gold  standard  is  normalized 
to  the  mean  square  error  due  to  motion  artifacts  from  the 
conventional  3D  acquisition. 

III.  RESULTS 
A.  Phantom  study 

7.  Trajectory  of  fiducial  in  projection  space  and 
phase  sorting 

In  Figs.  2(a)  and  2(b),  we  demonstrate  how  the  marker 
searching  algorithm  works.  The  intensity  of  a  projection  im¬ 
age  of  the  phantom  is  shown  in  Fig.  2(a)  for  the  selected 
ROI.  It  is  found  that  the  minimum  intensity  in  the  figure  did 
not  represent  the  marker.  After  filtering  and  subtraction,  the 
processed  ROI  image  is  shown  in  Fig.  2(b),  where  we  see 
that  the  local  minimum  gives  exactly  the  marker  location  on 
the  detection  plane.  The  fiducial  marker  was  identified  auto¬ 
matically  for  over  15  000  projections  accumulated  in  the 
phantom  experiments.  As  an  example,  the  coordinates  of  the 
marker  in  SI  direction  were  recorded  for  one  rotation  of  the 
10  mA  acquisition  in  the  first  experiment,  and  the  marker 
positions  in  the  lab  (real  world)  coordinate  system  were  cal¬ 
culated  as  well.  The  results  are  plotted  in  Fig.  3.  Although 
the  cone  beam  projected  the  marker  into  different  locations 
on  the  detection  plane  when  the  gantry  rotated  to  different 
angles,  it  is  found  that  the  peaks/valleys  of  the  projection 
positions  and  actual  positions  of  the  marker  correlated  well, 
which  suggests  the  marker  projection  positions  can  be  used 
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Fig.  2.  The  intensity  images  of  a  projection  containing  the  fiducial  marker 
before  (a)  and  after  (b)  processing.  After  the  image  processing,  the  minimum 
point  reflects  the  true  position  of  the  fiducial  marker  in  the  projection  space. 

to  determine  the  phase  directly;  the  amplitudes  in  the  projec¬ 
tion  images,  however,  do  not  reflect  the  real  displacements  of 
the  marker  and  should  not  be  used  for  the  data  binning. 

2.  Image  quality  as  a  function  of  number  of 
scans 

After  the  phase  of  each  projection  was  determined,  the 
projections  of  similar  phase  were  grouped  together  to  recon¬ 
struct  a  “single-phase”  image  using  the  conventional  FDK 
algorithm.  In  the  first  phantom  experiment,  we  studied  the 
mid-expiration  and  end-expiration  phases  with  6%  phase  un¬ 
certainty.  Results  from  mid-inspiration  and  end-inspiration 
phases  would  be  the  same  since  the  motion  was  sinusoidal 
(note  this  is  not  true  for  a  patient).  Since  the  phantom  had  a 
sinusoidal  motion  with  a  peak-valley  displacement  of 
2.0  cm,  the  mid-expiration  and  the  end-expiration  groups  of 
projections  had  residual  motion  of  about  0.375  and 
0.018  cm,  respectively.  Images  for  the  two  phases  were  re¬ 
constructed  for  a  series  of  numbers  of  scans,  and  are  shown 
in  Figs.  4  and  5.  To  save  space,  only  four  out  of  the  eight 
rotations  data  are  presented  here,  namely,  the  images  recon¬ 
structed  with  1 -rotations,  3-rotations,  5-rotations,  and 
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Fig.  3.  Automatically  extracted  respi¬ 
ration  signal  for  one  of  the  eight- 
rotations  in  the  10  mA  scan  of  the  mo¬ 
tion  phantom.  The  red  open  circles  are 
the  SI  coordinates  in  the  projection, 
and  the  solid  black  squares  are  the  SI 
coordinates  in  the  real  world  reference 
frame. 


7 -rotations  amount  of  projections.  Figure  6  shows  the 
“ground-truth”  images  at  the  two  phases.  As  expected,  the 
image  quality  gets  better  when  the  number  of  projections 
increases. 


The  relative  error  metric  is  plotted  in  Fig.  7  for  the  mid¬ 
expiration  and  end-expiration  phases.  Although  the  residual 
motion  amplitudes  in  the  two  phases  are  different,  the  rela¬ 
tions  between  the  image  quality  and  the  number  of  rotations 


Fig.  4.  Four-dimensional  single-phase  images  at  mid-expiration  phase  re¬ 
constructed  with  projection  data  exacted  from  one-turn,  three-turn,  five-turn, 
and  seven-turn,  respectively. 


Fig.  5.  Four-dimensional  single-phase  images  at  end-expiration  phase  re¬ 
constructed  with  projection  data  extracted  from  one-turn,  three-turn,  five- 
turn,  and  seven-turn,  respectively. 
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Fig.  6.  Stationary  reconstruction  at  mid-expiration  position  with  80  and 
10  mA  current  using  the  first  motion  mode. 


Fig.  7.  Relation  of  relative  error  and  the  number  of  rotations  for  mid¬ 
expiration  phase  (top)  and  end- expiration  phase  (bottom). 


Fig.  8.  Four-dimensional  CBCT  images  obtained  with  different  settings  re¬ 
sulting  in  the  same  radiation  dose.  From  top  to  bottom  row  are  the  images 
corresponding  to  80  mA  one  turn,  40  mA  two  turns,  20  mA  four  turns,  and 
10  mA  eight  turns,  respectively. 

have  the  same  trend.  The  reduction  of  the  relative  errors 
becomes  less  pronounced  as  the  number  of  rotations  in¬ 
creases.  After  4-5  rotations,  the  benefit  resulting  from  more 
rotations  starts  saturating. 

3 .  Tube  current  vs  number  of  scans 

In  the  second  experiment,  we  compare  images  obtained  at 
the  same  radiation  dose  level.  In  Fig.  8,  single-phase  images 
at  the  end-expiration  phase  are  shown  for  acquisition  settings 
of  one-rotation  80  mA,  two-rotation  40  mA,  four-rotation 
20  mA,  and  eight-rotation  10  mA.  Since  the  radiation  dose  is 
almost  linear  to  the  x-ray  tube  current  and  number  of  scans, 
these  acquisition  settings  have  basically  the  same  dose  level. 
The  images  of  3D  acquisitions  with  80  mA  for  the  phantom 
with  and  without  motion  are  shown  in  Fig.  9.  Again,  RE  is 
calculated  to  quantify  the  difference  of  the  images.  We  com¬ 
pare  the  REs  for  three  sets  of  single-phase  images,  each  hav¬ 
ing  a  different  phase  error  or  a  different  amount  of  residual 
motion.  For  example,  if  the  phase  interval  is  10%,  then  the 
image  phase  may  differ  from  the  prescribed  phase  by  5%, 
representing  the  phase  error  of  the  4D  image.  As  shown  in 
Fig.  10,  the  image  of  80  mA  acquisition  has  the  largest  rela¬ 
tive  error,  and  the  image  of  10  mA  has  the  smallest  relative 
error  value  for  all  three  curves.  For  the  same  patient  dose,  the 
low-current  and  multiple-rotation  acquisition  has  more  pro- 
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Fig.  1 1 .  Automatically  extracted  respiration  signal  of  the  lung  patient. 


Fig.  9.  Three-dimensional  CBCT  scan  of  the  phantom  with  (top)  and  with¬ 
out  (bottom)  motion  (second  mode)  with  80  mA  tube  current. 


jections,  and  thereby  more  data  can  be  generated  for  each 
sorted  single-phase  reconstruction.  From  the  curves  in  Fig. 
10,  it  is  seen  that  acquisition  with  more  rotations  and  lower 
currents  is  superior  to  high-current  acquisition  with  fewer 
rotations  in  terms  of  reconstruction  accuracy.  It  is  worth 
mentioning  that,  however,  the  metric  used  in  this  work  dif¬ 
fers  from  the  CNR  or  signal-to-noise  ratio  (SNR),  which 
may  behave  differently  with  the  tube  current. 

B.  Patient  study 

Figure  11  shows  the  extracted  respiration  signal  of  the 
lung  patient  from  his  CBCT  projections,  where  the  maxima 
represent  the  peak-inspiration  phase  and  the  minima  repre¬ 
sent  the  peak-expiration  phase.  As  mentioned  in  Sec.  II  B, 
the  extracted  signal  does  not  describe  true  displacements  of 
the  patient  anatomy;  however,  it  does  yield  the  correct  phase 
information. 


Fig.  10.  Image  error  of  4D  CBCT  at  different  settings  with  the  patient 
radiation  dose  kept  at  the  same  value.  The  image  phase  being  compared  is 
the  end  expiration.  The  phase  errors  are  5%  (dotted  line),  10%  (dot-dash 
line),  and  20%  (solid  line),  respectively.  A  larger  phase  uncertainty  means  a 
greater  residual  motion.  The  errors  are  normalized  to  that  of  the  3D  CBCT 
scan. 


Fig.  12.  Three-dimensional  and  4D  CBCT  images  of  the  lung  patient  shown 
in  transverse,  coronal,  and  sagittal  planes.  First  row  is  the  3D  CBCT  images, 
second  to  the  last  rows  are  end-inspiration  (phase  0%),  mid-expiration 
(phase  25%),  end- expiration  (phase  50%),  mid-inspiration  phase  (phase 
75%),  respectively. 
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In  Fig.  12,  we  show  the  reconstructed  images  in  trans¬ 
verse,  coronal,  and  sagittal  planes  for  the  3D  CBCT  data 
(one  of  the  four-turn  acquisition)  as  well  as  the  four  phases 
of  the  4D  data.  The  voxel  dimensions  in  the  reconstructed 
images  are  0.5  X  0.5  X  1.0  mm3.  It  is  observed  that  the  shape 
of  the  moving  structures  can  be  identified  more  accurately 
and  the  scan  provides  information  on  the  trajectory  of  these 
structures,  which  is  absent  in  the  3D  data.  Furthermore,  the 
view  aliasing  artifacts  appearing  in  Ref.  20  are  almost  elimi¬ 
nated  by  our  multiple-rotation  scheme. 

IV.  DISCUSSION  AND  CONCLUSION 

When  CBCT  is  used  in  imaging  the  thorax  or  upper  ab¬ 
domen  of  a  patient,  respiration  induced  artifacts  such  as  blur¬ 
ring,  doubling,  streaking,  and  distortion  are  often  observed, 
which  heavily  degrade  the  image  quality,  affecting  the  target 
localization  ability,  as  well  as  the  accuracy  of  irradiation 
dose  verification.  These  artifacts  are  caused  by  the 
respiration-induced  inconsistency  in  the  acquired  projection 
data.  Phase  binning  used  in  this  work  reduces  the  inconsis¬ 
tency  by  grouping  those  projections  according  to  their  phases 
and  removes/reduces  the  motion  artifacts.  To  obtain  high 
quality  4D  CBCT  images,  each  phase  must  contain  a  suffi¬ 
cient  number  of  projections.  It  is  known  that,  theoretically, 
the  spatial  resolution  that  can  be  resolved  in  a  reconstructed 
image  is  proportional  to  the  number  of  acquired  projections. 
More  precisely,  if  D  is  the  FOV  diameter,  r  is  the  resolution, 
then  N=ttDIv  is  the  number  of  projections  that  should  be 
provided  in  order  to  completely  restore  the  object  without 
noticeable  artifacts.  If  single-rotation  3D  CBCT  projections 
are  used  for  4D  phase-binned  CBCT  reconstruction,  dramatic 
artifacts  appear  because  of  the  limited  number  of 
projections.24  The  artifacts  can  be  clearly  seen  in  Figs.  4,  5, 
and  8.  We  have  shown  that  the  problem  can  be  solved  or  at 
least  greatly  relieved  by  multiple  gantry  rotations  with  the 
x-ray  tube  current  set  to  a  lower  value.  While  the  patient 
radiation  dose  remains  at  the  3D  CBCT  level,  the  new  4D 
CBCT  acquisition  protocol  provides  phase-resolved  images 
and  significantly  improves  the  quality  of  images. 

Radiation  dose  is  one  of  the  major  concerns  in  CT 
imaging,30  in  particular  the  on-board  CBCT  imaging  because 
of  the  repeated  use  of  modality  for  a  given  patient.  With 
increased  number  of  rotations  for  4D  acquisition,  we  found 
that  the  x-ray  tube  current  can  be  lowered  accordingly  with¬ 
out  compromising  the  image  quality.  In  a  sense,  this  radia¬ 
tion  dose  reduction  scheme  is  to  “spread”  the  photons  of  a 
projection  in  3D  CBCT  to  a  range  of  angles  for  4D  CBCT 
reconstruction,  which  represents  a  better  trade-off  between 
the  SNR  and  the  reduction  in  motion  artifacts.  In  the  pres¬ 
ence  of  significant  artifacts  (as  often  the  case  in  4D  CBCT 
imaging),  this  is  usually  preferable.  Of  course,  given  the 
same  number  of  rotations  (or  the  same  number  of  projections 
in  each  phase),  increasing  the  x-ray  tube  current  will  improve 
the  SNR  of  the  4D  images. 

In  addition  to  the  reduced  patient  dose,  the  heat  genera¬ 
tion  of  the  x-ray  tube  is  also  improved,  which  may  have 
some  practical  implication.  Indeed,  effective  removal  of  heat 


in  the  x-ray  tube  housing  has  been  problematic  in  on-board 
imager  based  CBCT  due  to  the  long  image  acquisition  time 
and  other  mechanical  difficulties.  In  the  current  implementa¬ 
tion  of  the  Trilogy™  system,  for  example,  the  lack  of  oil 
cooling  of  the  x-ray  tube  causes  rapid  increase  of  tempera¬ 
ture  when  the  beam  is  on  and  reaches  to  the  x-ray  housing 
temperature  threshold  after  three  continuous  scans. 

As  to  the  question  of  how  many  rotations  should  be  used 
for  the  4D  CBCT,  one  practical  consideration  is  the  total 
acquisition  time.  Although  theoretically  more  rotations  are 
preferable  because  it  gives  better  temporal  and  spatial  reso¬ 
lution,  it  prolongs  the  scanning  time  and  may  discomfort  the 
patient.  Possible  movement  of  the  patient  (other  than  the 
respiration)  may  cause  additional  image  artifacts  if  excessive 
number  of  rotations  is  used.  Our  experience  indicates  that 
acquisition  with  3-5  turns  at  10-30  mA  is  adequate  to  pro¬ 
vide  clinically  meaningful  4D  CBCT  images.  In  practice,  the 
decision  should  be  made  by  balancing  various  competing 
factors,  such  as  the  patient  radiation  dose,  scan  time,  and 
image  quality. 

An  alternative  strategy  to  increase  the  number  of  projec¬ 
tions  for  4D  CBCT  is  to  slow  down  the  gantry  rotation  dur¬ 
ing  the  data  acquisition.  It  has  the  advantage  of  less  user 
intervention  during  the  scanning  process,  since  all  required 
data  may  be  completed  in  a  single  rotation.  In  addition,  for 
multiple-rotation  4D  acquisition,  it  is  possible  that  there  is 
phase  overlapping  (projections  from  different  rotations  hav¬ 
ing  the  same  phase).  The  situation  can  be  relieved  to  a  cer¬ 
tain  extent  by  forcing  a  phase  shift  at  the  start  of  each  rota¬ 
tion.  In  clinic  practice,  for  example,  during  a  gated  radiation 
treatment,  the  patient  respiratory  motion  signal  can  be  moni¬ 
tored  through  a  real-time  tracking  system  (it  is  not  necessary 
for  the  system  to  be  synchronized  with  CBCT  in  our  work), 
and  the  acquisition  for  different  rotation  can  be  started  at 
different  phases  by  observing  the  motion  signal.  In  reality,  it 
was  also  found  that  the  phase  shift  is  not  critical,  primarily 
because  of  the  relatively  low  probability  for  all  scans  starting 
from  an  identical  phase  point  and,  more  or  less  irregular 
breathing  pattern  of  the  patients.  These  two  factors  reduce 
the  possibility  of  redundant  projections.  In  principle,  a  hy¬ 
brid  approach  combining  the  above  two  strategies  may  be  a 
viable  choice.  It  is  of  interest  and  significant  practical  impli¬ 
cation  to  systematically  study  the  different  data  acquisition 
scheme  and  find  the  optimal  protocol  for  clinical  application. 
This  work  entails  modification  of  the  manufacture  provided 
gantry  control  software  and  is  still  in  progress. 

In  summary,  we  have  demonstrated  that  the  commercially 
available  OBI  system  can  be  utilized  to  acquire  4D  CBCT 
images  with  reduced  motion/view- aliasing  artifacts  and  with¬ 
out  increasing  the  patient  radiation  dose.  To  achieve  this,  two 
steps  were  involved:  phase  binning  of  CBCT  projections  and 
CBCT  reconstruction.  The  phase  binning  in  our  approach 
was  accomplished  with  the  use  of  an  external  fiducial  marker 
and  the  analysis  of  its  trace  in  the  projection  space.  The 
binning  should  also  be  achievable  by  monitoring  an  im¬ 
planted  fiducial  or  even  an  anatomical  feature.  The  recon¬ 
struction  was  done  using  the  conventional  Feldkamp  algo¬ 
rithm  for  each  phase.  The  demand  for  more  projections  in 
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multi-phase  4D  CBCT  reconstruction  was  achieved  by  data 
acquisition  with  multiple  gantry  rotations.  We  found  that 
scaling  down  the  tube  current  in  accordance  with  the  in¬ 
creased  number  of  gantry  rotations  in  4D  CBCT  does  not 
compromise  the  quality  of  the  phased  images,  yet  has  the 
benefit  of  keeping  the  patient  radiation  dose  at  the  3D  CBCT 
level.  In  reality,  the  long  acquisition  time  of  this  4D  technol¬ 
ogy  may  have  limit  in  applications  such  as  real-time  thera¬ 
peutic  guidance.  However,  it  is  useful  in  patient  setup  and 
dose  verification  of  radiation  treatment  because  it  provides  a 
valuable  3D  or  4D  geometric  model  of  the  patient  just  prior 
to  treatment.  The  4D  CBCT  will  not  only  allow  us  to  better 
see  the  anatomy  by  eliminating  the  motion  artifacts,  but  also 
provides  kinetic  information  of  the  patient  anatomy,  which 
should  be  valuable  for  pre-treatment  verification  of  future  4D 
treatment  such  as  gated  radiotherapy  or  tumor  tracking. 
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Positron  emission  tonography  (PET)  is  useful  in  diagnosis  and  radiation  treatment  planning  for  a 
variety  of  cancers.  For  patients  with  cancers  in  thoracic  or  upper  abdominal  region,  the  respiratory 
motion  produces  large  distortions  in  the  tumor  shape  and  size,  affecting  the  accuracy  in  both 
diagnosis  and  treatment.  Four-dimensional  (4D)  (gated)  PET  aims  to  reduce  the  motion  artifacts 
and  to  provide  accurate  measurement  of  the  tumor  volume  and  the  tracer  concentration.  A  major 
issue  in  4D  PET  is  the  lack  of  statistics.  Since  the  collected  photons  are  divided  into  several  frames 
in  the  4D  PET  scan,  the  quality  of  each  reconstructed  frame  degrades  as  the  number  of  frames 
increases.  The  increased  noise  in  each  frame  heavily  degrades  the  quantitative  accuracy  of  the  PET 
imaging.  In  this  work,  we  propose  a  method  to  enhance  the  performance  of  4D  PET  by  developing 
a  new  technique  of  4D  PET  reconstruction  with  incorporation  of  an  organ  motion  model  derived 
from  4D-CT  images.  The  method  is  based  on  the  well-known  maximum-likelihood  expectation- 
maximization  (ML-EM)  algorithm.  During  the  processes  of  forward-  and  backward-projection  in 
the  ML-EM  iterations,  all  projection  data  acquired  at  different  phases  are  combined  together  to 
update  the  emission  map  with  the  aid  of  deformable  model,  the  statistics  is  therefore  greatly 
improved.  The  proposed  algorithm  was  first  evaluated  with  computer  simulations  using  a  math¬ 
ematical  dynamic  phantom.  Experiment  with  a  moving  physical  phantom  was  then  carried  out  to 
demonstrate  the  accuracy  of  the  proposed  method  and  the  increase  of  signal-to-noise  ratio  over 
three-dimensional  PET.  Finally,  the  4D  PET  reconstruction  was  applied  to  a  patient  case.  ©  2006 
American  Association  of  Physicists  in  Medicine.  [DOI:  10.1118/1.2192581] 


I.  INTRODUCTION 

Positron  emission  tomography  (PET)  scan  has  been  used  for 
effective  diagnosis,  staging,  and  evaluation  of  response  to 
chemo-  and  radiotherapy  for  a  variety  of  cancers,  such  as 
malignant  lymphoma,  malignant  melanoma,  lung  cancer, 
gastrointestinal  cancer,  and  head-and-neck  cancer.  For  lung 
cancer,  for  example,  by  combining  CT  information  with 
fluorodeoxyglucose(FDG)-PET,  it  helps  differentiate  tumors 
from  atelectasis  and  pleural  effusion,1  and  can  thus  improve 
the  delineation  of  the  tumor.  In  reality,  a  typical  PET  scan 
needs  3-7  min  at  each  field  of  view  (FOV),  the  acquired 
data  are  thus  averaged  over  many  breathing  cycles.  During 
the  long  acquisition  time,  the  lung  tumors  with  free  breathing 
can  move  up  to  30  mm,2  resulting  in  motion-blurring  of  the 
tumor  activity  in  the  PET  images  and  consequently  uncer¬ 
tainties  in  the  size,  position,  and  shape  of  the  tumor.3,4  It  has 
been  reported  that  the  respiratory  motion  can  lead  to  lower 
contrast  and  overestimation  of  the  lesion  volume.5  In  addi¬ 
tion,  the  respiration  motion  will  also  affect  the  co¬ 
registration  of  PET-CT  images  due  to  different  breathing 
conditions  of  the  two  modalities  during  data  acquisition.6 
Although  it  is  considered  as  one  of  the  important  criteria  in 
characterizing  the  lesions,  quantification  of  tracer  concentra¬ 
tion  cannot  be  achieved  accurately  when  conventional  PET  is 
used  for  thoracic  or  upper  abdominal  imaging. 

One  method  for  combating  the  respiration-induced  blur¬ 
ring  in  PET  is  to  perform  a  gated  four-dimensional  (4D) 
acquisition.  In  this  method,  the  PET  acquisition  is  triggered 
with  a  signal  characterizing  the  breathing  cycle  such  as 


RPM,  which  is  similar  to  the  gated  cardiac  PET  triggered  by 
the  ECG  signals.  The  acquired  data  are  divided  into  sepa¬ 
rate  storage  locations  based  on  the  respiratory  phase,  and 
cyclically  summed  over  many  respiratory  cycles.  Typically, 
the  4D  PET  acquisition  divides  each  respiratory  cycle  into 
0.2-1  s  time  bins  and  the  data  are  dumped  into  the  appro¬ 
priate  time  bin  during  the  entire  acquisition  period,  ranging 
from  5  to  10  min.  When  the  data  in  each  bin  are  recon¬ 
structed,  a  series  of  volumetric  PET  images  are  obtained, 
each  corresponding  to  a  respiration  phase  point.  4D  PET 
affords  useful  kinetic  information  of  organs  and  tumors,  and 
improves  the  quantification  of  tumor  volume  measurement.10 

While  the  gating  acquisition  scheme  may  reduce  tumor 
distortion  and  improve  tumor  volume  and  tracer  concentra¬ 
tion  measurements,  4D  PET  so  obtained  has  low  statistics 
and  generally  prevents  high-accuracy  quantification.  This  is 
because  a  similar  amount  of  projection  data  as  in  a  3D  PET 
scan  is  divided  into  many  different  time  frames,  and  the 
number  of  coincidence  events  per  image  frame  is  consider¬ 
ably  reduced,  leading  to  higher  image  noise  compared  to  the 
static  three-dimensional  (3D)  case.  A  choice  therefore  must 
be  made  between  short-duration  frames,  which  are  noisy  but 
preserve  temporal  resolution,  or  else  longer  frames  which  are 
statistically  superior  but  lose  temporal  resolution. 

Many  efforts  have  been  made  to  improve  the  statistics  of 
4D  PET  while  conserving  its  temporal  resolution.11-13  For 
example,  Rahmim  et  al.  proposed  to  incorporate  motion  in¬ 
formation  into  expectation-maximization  (EM)  reconstruc¬ 
tion  algorithm  to  compensate  for  the  rigid  whole-body 
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Fig.  1.  Flow  chart  of  the  4D  model- 
based  reconstruction  algorithm. 


motion.11  Qiao  et  al.  derived  from  first  principles  a  general 
list-mode  reconstruction  algorithm  that  accounts  for  nonrigid 
motion.  ’  For  cardiac  gated  PET,  Klein  et  al.  developed  a 
4D  deformable  registration  technique,  which  was  used  to 
sum  PET  images  of  different  frames  together  hence  reducing 
the  noise  in  the  resulted  image.14  Thorndyke  et  al.  proposed 
a  retrospective  stacking  technique  that  allows  one  to  enhance 
the  quality  of  respiration  gated  PET  images  by  registering 
different  breathing  phases  using  a  B Spline  model.15  Liviera- 
tos  et  al.  applied  rigid-body  transformation  to  list-mode  car¬ 
diac  PET  projection  data  to  correct  respiration  induced 
motion.16  In  this  work,  we  provide  an  alternative  method 
which  integrates  a  patient-specific  motion  model  into  PET 
reconstruction  process  and  develop  a  model-based  recon¬ 
struction  algorithm  to  achieve  a  better  noise  property  while 
conserving  the  temporal  resolution  for  4D  respiration-gated 
PET  imaging.  An  iterative  scheme  is  proposed  for  the  recon¬ 
struction  that  is  analogous  to  the  maximum-likelihood 
expectation-maximization  (ML-EM)  algorithm,17,18  in  which 


all  projection  data  acquired  at  different  phases  are  combined 
together  to  update  the  emission  map  with  the  use  of  a  respi¬ 
ration  motion  model.  The  motion  model  in  our  study  is  de¬ 
rived  from  4D  CT  images,  based  on  an  observation  that  4D 
PET  and  4D  CT  images  correlate  well  with  a  modem  com¬ 
bined  PET/CT  scanner.19  During  the  revision  of  this  paper, 
we  became  aware  that  a  similar  idea  has  been  presented  in 
Ref  13. 

The  proposed  method  is  validated  with  both  computer 
simulation  and  physical  phantom  studies.  In  Sec.  II,  we  first 
briefly  describe  the  4D  PET/CT  acquisition  procedure  and  a 
nonrigid  registration  technique  for  extracting  the  respiration 
motion  model  from  4D  CT  images,  then  the  model-based  4D 
reconstmction  algorithm  is  presented.  In  Sec.  Ill,  the  conver¬ 
gence  of  the  proposed  algorithm  is  studied,  and  the  recon¬ 
structed  result  is  compared  with  those  from  the  conventional 
3D  and  4D  methods.  Furthermore,  the  developed  algorithm 
is  applied  to  a  4D  PET  scan  of  a  physical  phantom.  We 
conclude  in  Sec.  IV. 
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Fig.  2.  Attenuation  and  emission  map 
for  the  PET  simulations  and  their  eight 
deformed  phases. 
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II.  METHOD 

A.  4D  PET/CT  acquisition 

4D  PET/CT  are  acquired  with  a  combined  PET/CT  scan¬ 
ner  (GE  Medical  Systems,  MI),  and  correlated  through  a  res¬ 
piratory  motion  tracking  system,  for  example,  the  real-time 

0.5  mCi  FDG 


1.5  cm 


Fig.  3.  Phantom  for  4D  PET/CT  experimental  study. 


position  management  (RPM)  system  (Varian  Medical  Sys¬ 
tems,  Palo  Alto,  CA),  which  interfaces  with  the  PET/CT 
scanner  so  that  the  acquired  data  are  related  to  the  real-time 
respiratory  motion  data.  For  4D  CT  acquisition,  one  scanning 
protocol  that  has  recently  been  developed  uses  a  “step-and- 
shoot”  technique  (cine  mode),  in  which  CT  projections  are 
acquired  repeatedly  over  a  complete  respiratory  cycle  at  each 
couch  position.  ’  The  period  of  each  CT  acquisition  seg¬ 
ment  (i.e.,  each  rotation)  is  time  stamped  with  an  “x-ray  ON” 
signal  and  recorded  by  the  position-monitoring  system.  Then 
4D  CT  images  are  retrospectively  sorted  into  groups  accord¬ 
ing  to  the  recorded  motion  data.  For  4D  PET  acquisition,  the 
photons  are  collected  prospectively  in  a  “gated  mode,”  simi¬ 
lar  to  the  cardiac  gating  technique.7-9  The  PET  system  syn¬ 
chronizes  data  acquisition  with  respiratory  motion  tracking 
system,  which  is  set  to  provide  triggers  (mimicking  the  ECG 
pulses  used  in  cardiac  gating)  at  the  end  inspiration  of  each 
breathing  cycle,  and  the  PET  system  starts  to  acquire  data 
into  several  0.2-0. 5  s  time  bins.  Total  acquisition  time  is 
typically  several  minutes.  Since  the  phase  sorting  techniques 
of  4D  CT  and  4D  PET  are  different  (retrospective  or  pro¬ 
spective),  to  mach  the  phase  of  4D  CT,  we  select  the  nearest 
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Fig.  4.  3D  PET  (top)  and  4D  PET  images  (middle  and  bottom)  at  phase  0. 
The  middle  is  the  regular  4D  PET,  and  the  bottom  is  image  from  4D 
deformation-based  reconstruction. 


bin  from  PET  data  by  analyzing  the  average  respiratory  pe¬ 
riod  according  to  its  RPM  respiration  trace,  more  details  on 

22  23 

this  method  can  be  found  in  Nehmeh  et  al.  ’ 

B.  Derivation  of  motion  model  from  4D-CT  images 

The  combined  PET/CT  scanner  provides  a  hardware  reg¬ 
istration  of  the  PET  and  CT  images  by  aligning  the  spatial 
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Fig.  5.  (a)  Convergence  of  the  4D  model-based  reconstruction  algorithm; 
(b)  comparison  of  the  total  uptake  value  in  the  tumor,  where  the  cross 
represents  the  value  measured  with  4D  deformation-based  reconstruction, 
the  open  circle  is  for  the  regular  4D  PET,  and  the  closed  square  is  for  the  3D 
PET. 


coordinates  of  the  two  sets  of  images.  At  the  same  time,  the 
4D  PET  and  CT  images  can  also  be  registered  temporally  by 
selecting  the  same  respiratory  phases  out  of  them.  In  order 
to  derive  a  motion  model  from  the  4D  CT  images  applicable 
for  4D  PET  image  reconstruction,  the  CT  images  are  down- 
sampled  from  512X512  array  to  256X256  array,  and  re¬ 
binned  into  3.27  mm  slice  thickness  to  match  that  of  PET 
images.  The  rebinning  process  also  ensured  that  the  center 
locations  of  the  CT  images  in  each  phase  coincide  with  the 
center  locations  of  the  PET  images.  The  CT  images  are  gen¬ 
erally  of  better  image  quality  with  detailed  anatomical  infor¬ 
mation.  Thus  the  respiration-induced  motion  is  represented 
more  accurately  by  CT  data  than  PET  data. 

The  motion  is  extracted  from  the  4D  CT  data  using  de¬ 
formable  registration  model  to  account  for  the  complex  or¬ 
gan  movement  and  shape  change.  We  adopted  the  free-form 
Spline  model  (B Spline)25,26  in  this  work  to  register  different 
phases  of  the  CT  images.  The  simplicity  and  yet  accuracy  of 
the  B  Spline  method  make  it  a  useful  tool  for  many  clinical 
applications.  In  this  model,  a  lattice  of  user-defined 


Medical  Physics,  Vol.  33,  No.  5,  May  2006 


1292 


Li  et  al 4D  PET  reconstruction 


1292 


Fig.  6.  3D  ungated  and  four  phases  of 
the  gated  PET  data  in  transverse  (left), 
coronal  (middle),  and  sagittal  (right) 
views.  From  top  to  bottom  are  3D  un¬ 
gated  PET,  phase  0%  (end- 
inspiration),  25%  (mid- expiration), 
50%  (end-expiration),  and  75%  (mid¬ 
inspiration),  respectively.  (Note  that 
no  tumor  is  shown  in  the  transverse 
view  of  phase  0%  because  they  had 
moved  out  of  view  of  this  slice.) 


nodes  is  overlaid  on  the  image.  Each  node  associates  a  de¬ 
formation  vector,  whose  components  are  determined  by  an 
optimization  procedure.  The  deformation  at  any  point  of  the 
image  is  calculated  by  spline  interpolation  of  the  adjacent 
nodes  values.  An  advantage  of  the  B Spline  model  is  that  the 
nodes  are  locally  controlled,  thus  the  displacement  of  an  in¬ 
terpolation  point  is  influenced  only  by  the  nearest  grid  points 
and  changing  a  lattice  node  only  affects  the  transformation 
regionally,  making  it  efficient  in  describing  local  deforma¬ 
tions.  Suitable  node  deformations  are  obtained  using  the 
gradient-based  L-BFGS  algorithm,  which  iteratively  varies 
the  deformations  until  the  registration  metric,  a  mathematical 
measure  of  similarity  between  images,  is  minimized.  The 
normalized  cross  correlation  (NCC)  was  used  as  the  metric 
in  this  study,  which  is  defined  as 


NCC  =  - 


N 

^i=ixjyj 


S;,  (*,)%>/ 


(i) 


where  x,y  are  two  images.  With  deformable  registration,  all 
phases  were  registered  to  the  end-inspiration  phase  (phase 


0%),  resulting  in  a  series  of  transformations  representing  a 
temporal  sequence  of  3D  deformation  fields,  or  in  other 
words,  a  4D  model  of  organ  motion. 


C.  4D  model-based  reconstruction  algorithm 

Let  y\  be  the  photon  numbers  collected  in  the  tth  time  bin 
(or  in  other  words,  phase  t)  at  detector  i  (or  line  of  response), 
and  we  assume  it  follows  a  Poisson  distribution: 


P(y‘, \n- 


e  yiy  ■ 


y}- 


(2) 


J'=E44  (3) 

j 

where  p(y\\y  J)  is  the  probability  of  measuring  y\  photons,  \j 
is  the  amount  of  radioactivity  in  pixel  j  of  the  emission  map 
A  at  the  phase  t ,  and  a\j  is  the  probability  that  the  radioac¬ 
tivity  in  j  contributes  to  the  count  level  in  detector  i.  The 
photon  counts  due  to  randoms  or  scatter  is  assumed  to  be 
precorrected,  and  is  not  considered  in  the  derivation  of  the 


Medical  Physics,  Vol.  33,  No.  5,  May  2006 


1293 


Li  et  al 4D  PET  reconstruction 


1293 


Fig.  7.  An  example  of  deformation  field,  which  was  derived  from  mapping 
CT  images  of  phase  25%  to  phase  0%. 


reconstruction  algorithm.  Let  Dfx  denote  the  transformation 
(motion  model)  derived  by  registering  two  4D  CT  images  at 
phase  tl  and  t2 ,  then  the  emission  map  at  phase  t  can  be 
derived  from  phase  0  as:  A^DqA0,  or 

=  =  (4) 

n 

where  is  the  element  of  the  deformation  matrix  which 
transform  the  image  at  phase  0  to  image  of  phase  t.  Similarly, 
the  “inverse”  deformation  is  defined  as 

(5) 

n 

The  likelihood  of  the  4D  PET  acquisition  is 

L(Y,A)  =  Up(ym-  (6) 

i,t 

The  maximum  likelihood  reconstruction  seeks  the  image 
A0  that  maximizes  the  above  objective  function  (6),  or 
equivalently, 

LL(7,  A)  =  2  2  v-  ln(>>  ■)->>■ 

t  i 

=  2  2  y‘  . 

(7) 


Fig.  8.  Reconstructed  image  of  3D  PET  (top),  regular  4D  PET  (middle),  and 
4D  model-based  reconstruction  (bottom). 


In  the  following,  by  analog  to  the  EM  (expectation- 
maximization)  algorithm  derived  for  3D  emission 
tomography,  ’  we  present  an  iterative  algorithm  for  the 
solution.  In  next  section,  we  show  that  the  likelihood  of  Eq. 
(7)  increases  monotonically  with  iterations  by  computer 
simulation  study, 


(X0)(t+D  = 

2j,Jdpj^iaij  Un 


y  dt,oy  j - h. - 

"T  mZauymfr 


(8) 


(Ow  =  2<A,°)w  (9) 

j 

where  i  e  [1  ,A^]  and  N  is  the  number  of  detector  bins; 
j ,m,n,p  g  [1  ,M]  and  M  is  the  number  of  pixels;  t  e  [1  ,T] 
and  T  is  the  number  of  phases;  and  k  is  the  iteration  number. 
The  implementation  of  Eq.  (8)  can  be  summarized  as  the 
flow  chart  in  Fig.  1.  Since  the  projection  data  in  every  time 
bin  are  used  together  to  update  the  emission  map,  the  statis¬ 
tics  is  expected  to  be  improved  by  the  proposed  algorithm. 
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Fig.  9.  Vertical  profiles  across  the  tumor  (indicated  by  the  lines  in  Fig.  8), 
where  the  green  dash-dotted  line  is  for  3D  PET,  the  blue  dash  line  is  for 
regular  4D  PET,  the  red  dotted  line  is  for  the  proposed  4D  model-based 
reconstruction,  and  the  black  solid  line  is  for  the  3D  PET  without  phantom 
motion. 


D.  Computer  simulation  and  physical  phantom 
studies 

In  the  simulation  study,  a  digital  phantom  of  attenuation 
map  was  obtained  from  CT  image  of  a  patient,  converted  to 
attenuation  coefficients  for  511  keV  y-ray  (Fig.  2). 32  The 
emission  phantom  was  an  artificial  image  obtained  by  seg¬ 
menting  out  the  inner  soft  tissue  from  the  CT  image  and 
adding  a  circular  “tumor,”  as  shown  in  the  same  figure.  The 
array  sizes  of  both  phantoms  were  256  X  256,  and  the  diam¬ 
eter  of  the  tumor  was  5  cm  (25  pixels).  The  activity  was 
arbitrarily  set  as  2.25  in  the  tumor  and  1.65  in  the  surround¬ 
ing  tissues.  Both  the  attenuation  and  emission  maps  were 
deformed  into  eight  phases  with  eight  transformation  matri¬ 
ces  D*tl  that  were  known  exactly.  Then  4D  PET  projections 
were  simulated  with  the  deformed  maps  by  line  integrals. 
The  array  size  of  the  simulated  sinograms  was  256X256. 
Poisson  noise  was  subsequently  added  into  each  of  the  eight 
simulated  4D  PET  data  sets  using  random  number  generator, 
corresponding  to  a  level  of  10  000  photon  counts  per  view. 

A  physical  phantom  as  shown  in  Fig.  3  was  placed  on  the 
top  of  a  platform  capable  of  sinusoidal  motion  along  the 
cranial-caudal  direction.  The  amplitude  of  motion  was 
1.5  cm  with  a  period  of  4.3  s.  0.5  mCi  2-[18F]  fluoro-2- 
deoxy-D-glucose  (FDG)  isotope  was  injected  into  two  tubes 
inside  the  phantom,  each  with  size  of  0.75  cm2X3  cm,  re¬ 


sembling  two  tumors.  Another  5  mCi  isotope  was  injected 
into  the  phantom  as  the  background  activities,  and  diluted 
with  approximately  5000  ml  water.  A  clinical  4D  PET/CT 
protocol  was  used  to  scan  the  phantom.  Both  the  4D  CT  and 
4D  PET  covered  the  same  FOV  of  150  mm  in  axial  direc¬ 
tion.  A  Total  of  eight  phases  was  obtained  for  the  CT  images, 
and  30  0.2  s  bin  groups  were  acquired  for  the  4D  respiration¬ 
gated  PET  acquisition,  and  total  acquisition  time  of  the  4D 
PET  was  7  min. 

E.  Patient  study 

In  an  IRB -approved  study,  a  clinical  4D  PET/CT  scan  was 
performed  for  a  pancreatic  cancer  patient.  The  patient  re¬ 
ceived  16.5  mCi  18FDG  3  h  prior  to  the  scan,  whose  weight 
was  83  kg.  The  4D  CT  and  4D  PET  covered  the  same  FOV 
of  the  patient,  which  was  700X700  mm2  in  tranaxial  plane 
and  150  mm  in  axial  direction,  total  of  10  phases  was  ob¬ 
tained  for  the  CT  images,  and  again  30  0.2  s  bin  groups  were 
used  for  the  4D  respiration-gated  PET  acquisition.  The  maxi¬ 
mum  tumor  movement  in  the  patient  was  about  1.3  cm,  mea¬ 
sured  from  his  4D  CT  images.  The  acquisition  time  of  the  4D 
gated  PET  was  10  min. 

III.  RESULTS 

A.  Computer  simulation 

The  simulated  4D  PET  data  were  reconstructed  with  the 
proposed  method.  In  Fig.  4,  it  shows  the  reconstructed  image 
of  the  first  phase  with  the  proposed  method.  For  comparison, 
the  regular  4D  PET  reconstruction  result  using  one  single 
projection  data  set  was  also  presented,  as  well  as  the  image 
of  3D  PET,  which  was  obtained  by  reconstructing  the  aver¬ 
aged  projections  over  all  phases.  The  color  bar  associated 
with  each  picture  shows  the  intensity  of  the  presented  im¬ 
ages,  where  the  bright  color  indicates  higher  intensity  and 
the  dark  color  indicates  lower  intensity  with  arbitrary  units. 
The  algorithm  to  reconstruct  the  3D  PET  or  regular  4D  PET 
was  the  conventional  ML-EM  algorithm.  The  number  of  it¬ 
erations  was  40  for  all  three  methods.  From  these  pictures,  it 
is  easy  to  see  the  distortion  and  blurring  of  the  tumor  in  the 
3D-PET  image,  large  noise  in  the  single  phase  image  (regu¬ 
lar  4D  PET  image),  and  correctly  phase-resolved  image  with 
better  noise  control  using  the  proposed  4D  reconstruction 
technique. 


Table  I.  SNR  comparisons  for  regular  3D  PET  with  and  without  motion,  regular  4D  PET,  and  proposed  4D 
model-based  reconstruction. 


3D  PET 

with  motion 

4D  gated  PET 

4D  model-based 

reconstruction 

3D  PET 

without  motion 

Average  activities 

4.538±0.994 

8.377±2.468 

7.693  ±1.259 

7.985±0.984 

in  tumors 

Average  activities 

0.122±0.074 

0.179±0.134 

0.120±0.086 

0.119±0.076 

in  normal  tissues 

SNR 

4.43 

3.32 

6.00 

7.97 
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Fig.  10.  Examples  of  matched  phases 
of  PET  sinogram  data  (left  column) 
and  CT  (right  column).  From  top  to 
bottom  are  phase  0%,  10%,  20%,  and 
30%,  respectively,  with  0%  corre¬ 
sponding  to  the  end-inspiration  phase. 


To  examine  the  convergence  property  of  the  algorithm, 
the  objective  function  value  LL  of  Eq.  (6)  was  calculated  at 
each  iteration  and  plotted  in  Fig.  5(a)  for  100  iterations.  It  is 
found  that  the  objective  function  LL  increases  monotonically 
after  each  step  of  iterations,  and  similar  to  the  conventional 
EM  algorithm,  the  convergence  speed  gets  slower  after  a 
certain  number  of  iterations.  We  further  quantified  the  differ¬ 
ence  among  the  reconstruction  results  by  calculating  the  total 
uptake  value  (TUV),  which  was  defined  as  the  summation  of 
all  the  activities  in  the  tumor.  Since  the  tumor  region  was 


known  exactly,  quantitative  measurements  of  the  TUV  are 
equivalent  to  the  measurements  of  the  density  of  the  activity 
distribution  inside  the  tumor.  When  tumor  motion  is  present, 
the  radioactivities  span  a  larger  region  compared  to  the  static 
tumor.  Therefore,  3D  PET  acquisition  will  result  in  blurring 
of  the  tumor  boundary  and  decrease  of  the  TUV.  From  an¬ 
other  point  of  view,  accurate  TUV  indicates  accurate  mea¬ 
surement  of  the  tumor  size,  and  hence  reflects  the  preserving 
of  the  spatial  resolution  in  certain  degree.  In  Fig.  5(b),  the 
TUVs  were  plotted  versus  iteration  number  for  different 
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Fig.  11.  Reconstructed  images  of  3D 
PET  (top),  regular  4D  PET  (middle), 
and  4D  deformation-based  reconstruc¬ 
tion  for  the  pancreas  patient. 
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methods.  The  true  value  of  TUV  in  the  digital  phantom  was 
1107.3  with  arbitrary  unit.  As  expected,  the  3D  PET  under¬ 
estimated  the  TUV,  and  both  the  regular  4D  PET  and  the 
proposed  4D  method  give  a  similarly  better  estimate  of  the 
uptake  value,  which  shows  that  the  proposed  4D  model- 
based  reconstruction  method  preserves  the  quantitative  accu¬ 
racy  while  reduces  the  image  noise. 

B.  Phantom  experiment 

The  proposed  4D  reconstruction  algorithm  was  tested 
with  four  groups  (phases)  of  the  acquired  PET  data,  corre¬ 
sponding  to  end-inspiration  phase,  mid-expiration  phase, 
end-expiration  phase,  and  mid-inspiration  phase.  Their  coun¬ 
terpart  phases  of  CT  images  were  selected  for  respiratory 
motion  model  derivation.  Figure  6  shows  the  four  phases  of 
the  regular  4D  gated  PET  images  and  the  3D  ungated  PET  as 
well  in  the  transverse,  coronal,  and  sagittal  views.  The  defor¬ 
mation  fields  derived  from  the  4D  CT  images  by  B  Spline 
deformable  model  are  illustrated  in  Fig.  7  for  motion  of  end- 
inspiration  phase  to  mid-expiration  phase,  in  which  the  ar¬ 
rows  show  the  direction  and  amplitude  of  the  movement  of 
each  B Spline  node.  The  motion  derived  from  CT  image  reg¬ 
istration  was  close  to  the  real  phantom  movement,  however, 
noticeable  errors  still  existed. 

The  reconstruction  was  performed  using  the  proposed  al¬ 
gorithm  for  improved  4D  PET  imaging,  and  conventional 


ML-EM  algorithm  for  3D  PET  (with  and  without  motion) 
and  individual  phase  reconstruction  of  the  regular  4D  PET, 
as  done  in  the  computer  simulation  study.  Results  after  40 
iterations  are  shown  in  Fig.  8  for  one  sagittal  view.  The  pro¬ 
files  along  the  axial  direction  are  also  presented  in  Fig.  9  for 
one  of  the  tumors,  where  the  horizontal  axis  is  the  slice  num¬ 
ber  and  the  vertical  axis  is  the  relative  intensity  of  the  images 
with  arbitrary  units.  It  can  be  seen  that  the  regular  4D  gated 
PET,  the  3D  PET  without  motion  (which  is  our  “ground 
truth”),  and  the  model-based  4D  PET  have  similar  intensity 
profiles  and  tumor  locations  as  well,  indicating  that  our  pro¬ 
posed  algorithm  does  not  introduce  obvious  bias  into  the 
image  reconstruction  while  reducing  the  statistical  noise. 
From  both  Figs.  8  and  9,  we  can  find  that  cranial-caudal 
motion  resulted  in  the  broadening  of  the  tumor  in  3D  PET 
images:  the  tumor  volume  became  larger,  and  the  activity 
concentration  became  lower.  The  4D  respiration-gated  acqui¬ 
sition  can  successfully  suppress  the  tumor  blurring  and  re¬ 
store  the  tumor  uptake  value.  In  addition,  the  proposed  4D 
model-based  reconstruction  algorithm  preserved  the  tumor 
shape  compared  to  the  regular  4D  PET,  while  reducing  the 
image  noise  (see  Fig.  8).  Similar  results  can  be  also  found  for 
the  other  tumor. 

To  quantify  the  improvement,  we  further  measured  the 
activities  in  the  tumors  and  normal  tissues,  and  calculated  the 
signal-to-noise  ratio,  defined  as  SNR=\St-Sn\/y]o^+o^, 
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where  St  and  Sn  are  the  measured  average  activities  in  se¬ 
lected  regions  of  interest  (ROI)  in  the  tumors  and  normal 
tissues,  respectively,  and  the  crt  and  crn  are  their  standard 
deviations.  The  size  of  the  selected  ROI  in  this  study  was  96 
voxels  drawn  on  6  slices  for  both  tumor  and  normal  tissues, 
and  the  results  are  listed  in  Table  I  with  arbitrary  units.  Al¬ 
though  the  regular  4D  PET  presented  higher  signal  St,  it  is 
found  that  the  SNR  is  even  lower  than  3D  PET  due  to  the 
presence  of  strong  noise.  On  the  contrary,  the  proposed  4D 
reconstruction  method  increased  the  SNR  about  35%  over 
the  3D  PET. 

C.  Patient  study 

In  Fig.  10,  we  show  the  pre-corrected  PET  sinograms  to¬ 
gether  with  the  corresponding  CT  images  at  four  different 
phases,  namely,  phase  0%,  10%,  20%,  30%  (with  0%  corre¬ 
sponding  to  end-inspiration  phase).  The  reconstructed  PET 
images  using  different  methods  are  shown  in  Fig.  11,  where 
the  top  row  is  the  3D  ungated  PET  (obtained  by  putting  all 
the  acquired  4D-PET  projections  in  group),  the  second  row  is 
the  regular  4D  gated  PET,  and  the  last  row  is  the  model- 
based  4D  PET.  The  number  of  iterations  was  40  for  all  three 
reconstructions.  From  the  figure,  we  can  see  that  the  tumor  is 
better  defined  using  the  proposed  reconstruction  algorithm 
compared  to  the  regular  4D  gated  PET  due  to  less  contami¬ 
nation  of  the  statistic  noise.  The  SNRs  were  calculated  for 
the  three  images,  which  are  2.21,  1.83,  and  4.17  for  the  3D 
ungated  PET,  regular  4D  gated  PET,  and  model-based  4D 
PET,  respectively.  In  the  calculation,  the  signal  St  was  the 
mean  activity  of  the  central  125  voxels  of  the  pancreas  tu¬ 
mor,  and  the  background  noise  Sn  was  calculated  using  an 
arbitrarily  selected  uniform  region  containing  also  125  vox¬ 
els.  Note  that  the  phase-resolved  PET  sinogram  data  ob¬ 
tained  from  gated  acquisition  mode  is  actually  an  average  of 
similar  phases  from  different  respiration  cycles  (typically 
over  100  cycles),  it  is  expected  that  the  breath  motion  would 
not  have  the  perfect  repetition  as  of  the  phantom  study,  there¬ 
fore,  the  motion  model  derived  from  the  CT  phases  of  the 
patient  was  only  an  approximate  model  when  it  is  applied  to 
the  4D  gated  PET.  Nevertheless,  the  patient  study  shows  that 
the  improvement  of  the  image  quality  is  still  noticeable  by 
the  proposed  method. 

IV.  DISCUSSION  AND  CONCLUSION 

In  this  work,  we  have  demonstrated  the  feasibility  of  in¬ 
corporating  4D  motion  model  into  PET  image  reconstruc¬ 
tion,  and  our  computer  simulations  showed  the  property  and 
quantitative  accuracy  of  the  proposed  algorithm  given  the 
motion  model  is  exactly  known.  Although  the  physical  phan¬ 
tom  study  conducted  in  this  work  did  not  represent  the  true 
clinical  situation  (for  example,  it  was  only  a  rigid  one¬ 
dimensional  “whole-body”  motion),  it  is  useful  for  the  ro¬ 
bustness  test  of  the  proposed  method  as  well  as  the  quanti¬ 
tative  evaluation  of  the  accuracy.  The  motion  model  used  in 
this  work  is  derived  from  the  4D  CT  images  of  the  same 
subject.  Comparing  with  direct  registration  of  4D  PET 
frames,14  motion  model  derived  from  CT  frames  is  usually 


more  accurate  because  of  the  higher  spatial  resolution  and 
less  noise  in  CT  image.  An  alternative  approach  to  improving 
the  4D  PET  statistics  is  to  reconstruct  the  low-count  PET 
frames  first,  then  apply  the  CT-based  motion  model  for  ret¬ 
rospectively  stacking  similar  to  Ref.  14.  Systematic  compari¬ 
sons  among  these  alternative  methods  would  be  of  great  re¬ 
search  interest. 

Incorporating  the  motion  model  into  the  inverse  process  is 
a  general  approach  in  estimation  and  have  been  applied  to 
CT  and  MRI  reconstruction  ’  and  radiation  dose 
optimization.35  In  this  process,  the  motion  model  serves  as 
a  priori  knowledge  of  the  system  and  provides  valuable  par¬ 
tial  guidance  to  searching  the  solution.  The  method  can  be 
also  applied  to  4D  PET  acquired  with  list  mode.  A  potential 
shortcoming  of  the  proposed  reconstruction  method  (or  other 
model-based  stacking  methods)  is  that  it  does  not  consider 
the  possible  change  in  patient’s  organ  motion  pattern.  Indeed, 
the  motion  may  change  to  some  degree  from  the  time  of  4D 
CT  acquisition  to  the  gated  PET  acquisition,  which  is  done 
sequentially  instead  of  simultaneously  on  a  PET/CT  scanner. 
In  particular,  the  frame  obtained  in  4D  PET  is  actually  an 
average  of  the  similar  phases  of  different  respiration  cycles. 
Therefore,  the  model  derived  from  4D  CT  does  not  match 
exactly  the  PET  data.  However,  it  should  be  a  reasonably 
good  estimation.  Unless  the  organ  motion  pattern  is  com¬ 
pletely  reversed,  we  anticipate  that  the  proposed  technique 
will  help  to  reduce  the  motion  artifacts  even  if  a  slight 
change  in  the  motion  pattern  occurs.  In  practice,  a  consistent 
respiratory  pattern  should  be  maintained  through,  for  ex¬ 
ample,  a  proper  voice  coaching,  during  the  PET  data  acqui¬ 
sition  process.  This  research  also  highlights  the  importance 
of  the  development  of  a  robust  organ  motion  tracking  system 
in  the  future.  With  this,  one  may  obtain  the  4D  patient  model 
during  the  4D  PET  process  and  then  incorporate  it  into  the 
PET  image  reconstruction  process.  In  radiotherapy  clinics, 
4D  CT  has  been  adopted  rapidly  to  identify  the  tumor  motion 
and  to  guide  the  4D  treatment  planning  process.36  Using  the 
4D  CT  data  for  the  removal  of  motion  artifacts  in  PET  ac¬ 
quisition  seems  to  be  logical  and  fit  well  with  the  data  flow 
of  modern  image-guide  radiotherapy. 

In  summary,  we  have  proposed  a  new  reconstruction  al¬ 
gorithm  for  improved  4D  PET  imaging.  It  integrates  all  the 
projection  data  acquired  at  different  time  bins  into  one  ob¬ 
jective  function,  therefore,  avoiding  the  drawback  of  low- 
count  nature  of  regular  4D  gated  PET  acquisitions.  The  pro¬ 
posed  algorithm  seems  to  converge  monotonically  to  the 
solution  in  both  cases.  The  improvement  of  the  image  quality 
was  clearly  demonstrated.  In  the  simulation  study,  the  pro¬ 
posed  method  resulted  in  less  noise  in  the  reconstructed  im¬ 
age  while  keeping  the  accuracy  of  tumor  uptake  value  mea¬ 
surement.  When  applied  to  the  phantom  experiment,  the 
proposed  algorithm  led  to  an  increase  of  SNR  over  both 
regular  3D  and  4D  PET,  and  reduced  the  motion  artifacts. 
This  work  should  be  useful  in  improving  delineation  of  mov¬ 
ing  tumors  and  quantitative  accuracy  of  the  tumor  uptake 
measurements. 
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Abstract 

On-board  imager  (OBI)  based  cone-beam  computed  tomography  (CBCT)  has 
become  available  in  radiotherapy  clinics  to  accurately  identify  the  target  in 
the  treatment  position.  However,  due  to  the  relatively  slow  gantry  rotation 
(typically  about  60  s  for  a  full  360°  scan)  in  acquiring  the  CBCT  projection 
data,  the  patient’s  respiratory  motion  causes  serious  problems  such  as  blurring, 
doubling,  streaking  and  distortion  in  the  reconstructed  images,  which  heavily 
degrade  the  image  quality  and  the  target  localization.  In  this  work,  we  present  a 
motion  compensation  method  for  slow-rotating  CBCT  scans  by  incorporating 
into  image  reconstruction  a  patient-specific  motion  model,  which  is  derived 
from  previously  obtained  four-dimensional  (4D)  treatment  planning  CT  images 
of  the  same  patient  via  deformable  registration.  The  registration  of  the  4D 
CT  phases  results  in  transformations  representing  a  temporal  sequence  of 
three-dimensional  (3D)  deformation  fields,  or  in  other  words,  a  4D  model  of 
organ  motion.  The  algorithm  was  developed  heuristically  in  two-dimensional 
(2D)  parallel-beam  geometry  and  extended  to  3D  cone-beam  geometry.  By 
simulations  with  digital  phantoms  capable  of  translational  motion  and  other 
complex  motion,  we  demonstrated  that  the  algorithm  can  reduce  the  motion 
artefacts  locally,  and  restore  the  tumour  size  and  shape,  which  may  thereby 
improve  the  accuracy  of  target  localization  and  patient  positioning  when  CBCT 
is  used  as  the  treatment  guidance. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

A  new  technology  of  cone-beam  computed  tomography  (CBCT)  has  recently  been  integrated 
on-board  with  the  linear  accelerator  (linac)  in  radiotherapy  clinics.  Superior  to  the  common 
approaches  based  on  two  orthogonal  images  provided  by  the  mega-voltage  electronic 


0031-9155/06/020253+15$30.00  ©  2006  IOP  Publishing  Ltd  Printed  in  the  UK 


253 


254 


T  Li  et  al 


portal  imaging  device  (EPID),  CBCT  can  provide  high-resolution  three-dimensional  (3D) 
information  of  the  patient  anatomy  in  the  treatment  position,  and  thus  has  great  potential  for 
improved  target  localization  and  irradiation  dose  verification  in  radiotherapy  (Jaffray  et  al 
1999,  Moseley  et  al  2004,  Sidhu  et  al  2003,  van  Herk  et  al  2004),  and  can  also  be  utilized  in 
synchronized  respiratory  gating  radiotherapy  (Jin  and  Yin  2005).  However,  when  the  on-board 
CBCT  is  used  in  imaging  the  thorax  or  abdomen  of  a  patient,  respiration  induced  artefacts  such 
as  blurring,  doubling,  streaking  and  distortion  are  observed,  which  heavily  degrade  the  image 
quality,  and  affect  the  target  localization  ability,  as  well  as  the  accuracy  of  dose  verification 
(Sonke  et  al  2005).  These  artefacts  are  much  more  severe  than  those  found  in  conventional 
CT  examinations.  In  conventional  CT,  each  rotation  of  the  scan  can  be  completed  within 
1  s,  during  this  period  the  organ/tumour  motion  is  relatively  small.  Furthermore,  patient 
body-restraints  and  breath-hold  techniques  can  be  used  to  minimize  the  motion  if  necessary. 
In  contrast,  in  a  CBCT  scan,  the  gantry  rotation  speed  is  much  slower,  typically  40  s  to  1  min 
for  a  full  360°  scan  in  acquiring  the  projection  data,  which  covers  more  than  10  breathing 
cycles  for  most  patients.  Breath  holding  is  uncomfortable  or  even  impossible  for  someone 
such  as  paediatric  or  lung  cancer  patients.  The  large  and  complex  movement  of  organs  during 
the  data  acquisition  causes  much  more  serious  problems  in  CBCT  than  the  conventional  CT. 

Considerable  efforts  have  been  made  to  investigate  efficient  methods  to  reduce  the  motion 
artefacts  in  conventional  CT  and  other  imaging  modalities  such  as  magnetic  resonance  imaging 
(MRI)  and  positron  emission  tomography  (PET)  (Atalar  and  Onural  1991,  Buhler  et  al  2004, 
Crawford  et  al  1996,  Cuppen  et  al  1985,  Dhanantwari  et  al  2001,  Ritchie  et  al  1996,  Wang  and 
Vannier  1995,  Willis  and  Bresler  1995).  Wang  and  Vannier  (1995)  developed  a  patient  motion 
estimation  and  compensation  technique  for  helical  CT  systems,  and  showed  good  simulation 
results,  but  it  was  limited  to  translational  motion  of  the  whole  patient  and  did  not  extend  to 
organ  motion.  Willis  and  Bresler  (1995)  cast  the  motion  artefact  problem  as  a  time-varying 
tomography  problem  and  proposed  special-purpose  hardware  to  optimally  sample  the  spatially 
and  temporally  band-limited  CT  signal  space.  A  parametric  model  for  the  respiratory  motion 
was  first  used  in  MRI,  and  the  motion  artefacts  were  successfully  reduced  by  modifying 
the  reconstruction  algorithm  (Atalar  and  Onural  1991,  Cuppen  et  al  1985).  Crawford  et  al 
(1996)  brought  the  concept  into  CT  imaging  and  derived  an  exact  reconstruction  formula  for 
motion  compensation  for  CT  scans  for  parallel-beam  projections.  The  method  is  referred  to 
as  the  CTX  method.  In  this  approach,  the  respiratory  motion  was  based  on  magnification  and 
displacement  of  the  object,  and  the  backprojection  was  performed  in  a  reference  frame  that 
moved  according  to  the  motion  model.  Ritchie  et  al  (1996)  pointed  out  that  the  usefulness 
of  CTX  was  limited  by  the  fact  that  the  time-varying  magnification  model  was  not  valid  for 
motion  in  the  chest.  They  extended  the  method  with  a  physiologically  more  correct  model 
for  respiratory  motion  and  applied  the  CTX  time- varying  model  on  a  local  basis  using  pixel- 
specific  backprojection  (PSBP).  A  method  for  estimating  the  in-plane  motion  of  every  pixel 
in  the  image  at  the  time  each  projection  is  acquired  was  also  developed.  However,  the  method 
depended  heavily  on  estimating  the  in-plane  motion  at  individual  manually  determined  node 
points.  Motion  correction  algorithms  that  assume  a  motion  model  work  well  when  the  motion 
conforms  to  the  model  but  have  limited  success  when  it  does  not  (Ablitt  et  al  2004,  Crawford 
et  al  1996,  Cuppen  et  al  1985,  Linney  and  Gregson  2001).  Furthermore,  the  CTX-based 
methods  cannot  be  applied  to  deal  with  the  general  3D  motion  in  CT  imaging. 

Recently,  four-dimensional  (4D)  CT  has  gained  popularity  in  guiding  radiation  treatment 
in  order  to  explicitly  account  for  the  respiratory  motion  (Low  et  al  2003,  Pan  et  al  2004,  Rietzel 
et  al  2005,  Vedam  et  al  2003).  With  multiple  scans  at  each  patient  couch  position,  it  generates 
a  series  of  phase  images  with  respect  to  the  motion  data  acquired  by  a  real-time  positioning 
system  during  the  scan.  The  phases  can  be  used  to  derive  a  patient-specific  deformation  field , 
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which  accurately  models  the  motion  pattern  of  the  patient  (Brock  et  al  2003,  Schreibmann 
et  al  2005).  Efforts  are  being  made  to  use  the  4D  patient  model  for  time-resolved  radiation 
therapy  planning  (Keall  et  al  2005,  Trofimov  et  al  2005,  Webb  2005).  In  this  work,  we  explore 
the  feasibility  of  incorporating  the  4D  patient  data  into  the  image  reconstruction  process  to 
obtain  phase-resolved  CBCT  images.  Incorporation  of  a  customized  object  motion  model  in 
CT  image  reconstruction  has  been  investigated  previously  (Crawford  et  al  1996,  Ritchie  et  al 
1996)  with  the  goal  of  reducing  motion  artefacts.  With  the  use  of  the  motion  information 
derived  from  the  same  patient,  the  proposed  method  should  lead  to  a  more  accurate  and 
robust  solution  to  the  problem.  In  the  following,  a  modified  filtered  backprojection  (FBP) 
algorithm  is  developed  for  the  simplest  two-dimensional  (2D)  parallel-beam  geometry,  and 
validated  with  translational  and  more  complex  motion,  with  and  without  Gaussian  noise.  A 
modified  Feldkamp  algorithm  is  then  implemented  for  CBCT,  and  tested  with  a  more  realistic 
deformable  phantom  constructed  from  a  patient’s  4D  CT  images. 

2.  Methods  and  materials 

2.7.  Reconstruction  with  deformation  field 

For  2D  parallel  geometry,  the  projection  of  an  object  g(x,  y),  R(0,  p),  at  gantry  rotational 
position  0  and  projection  distance  p  is  given  by 

/oo  poo 

/  g(x,  y)8(x  cosO  +  y  sin#  —  p)  dx  dy,  (1) 

-oo  7  — oo 

where  <$(.)  is  the  Dirac  delta-function.  The  image  reconstruction  gives  a  band-limited 
estimation  of  the  object,  g#(x,  y)  as  follows, 

Y  p^TT  pOO 

gB(x,y)  =  -  /  R(9,  p)f(xcos0  +ysind  -  p)dpd0,  (2) 

Z  7 0  7-00 

where  f{s)  =  f°^o)  \co\ W(co)  exp(2TC)cos)  dco  is  the  filter  function.  Different  types  of  filters 
can  be  obtained  by  different  window  function  W {(d)  designs.  In  practice,  equation  (2)  is 
usually  implemented  by  the  following  ‘filtered  backprojection’  steps, 

/oo 

R(9,  p)f  (p'  —  p)  dp,  (3) 

-oo 

1  f2n  - 

=  -  /  R(0,  x  cos  6  +  y  sinO)  dO,  (4) 

2  Jo 

where  (3)  is  the  filtered  step  implemented  by  convolution  (or  can  be  alternatively  implemented 
by  Fourier  transform),  and  (4)  is  the  backprojection  step. 

Now  if  the  object  moves  during  the  scan,  the  projection  data  acquired  at  each  angle 
9 1 ,  i  =  0,  1,  2,  . . . ,  N  —  1,  actually  correspond  to  a  series  of  ‘different  objects’  got  (or  more 
precisely,  different  status  of  the  object),  which  can  be  described  as  the  object  at  the  first 
phase  go  (or  equivalently,  at  the  first  projection  angle)  being  deformed  by  a  time-dependent 
transformation  Tq.  (i.e.,  a  4D  motion  model),  so  that 

get(x,  y)  =  Te,[go](x,  y).  (5) 

We  assume  that  the  transformation  T <9.  is  known  for  each  projection  angle,  so  is  the  ‘inverse’ 
transformation  V0j,  go(x,  y)  =  [g0.](x,  y).  It  will  be  discussed  later  on  how  to  use 
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In  order  to  apply  the  motion  model  to  reconstruction,  we  rewrite  the  reconstruction 
formula  (4)  in  the  following  discrete  form, 

N 

8b(x,  y)  =  (x,  y),  (6) 

i  =  l 

where  got  (x,  y)  =  R(ft,xcosft  +  y  sin  ft)  is  essentially  the  backprojection  of  filtered 
projection  at  angle  ft,  and  can  be  regarded  as  an  intermediately  reconstructed  object  with 
a  single-angle  projection.  By  summing  all  these  intermediate  objects  together  from  all  angles, 
some  pixels  are  enhanced  and  some  pixels  are  cancelled  out,  and  a  reconstructed  image  can 
be  obtained.  However,  when  motion  is  present,  the  corresponding  pixels  of  the  intermediate 
objects  are  misplaced,  and  the  summation  will  result  in  blurring,  doubling  or  other  distortions. 
Similar  to  the  assumption  used  by  Ritchie  et  al  (1996),  one  can  assume  that  local  correction 
was  a  valid  approximation  in  CT  reconstruction  and  the  backprojection  can  be  performed  in  a 
deformed  reference.  Therefore,  we  propose  heuristically  to  deform  the  intermediate  objects 
to  the  same  phase  before  doing  the  summation: 

N 

gB(x,  v)  ~  ^  X)T g.lge^x,  y)].  (7) 

i= 1 

Note  that  got  (x,  y)  of  similar  phases  can  be  grouped  together  before  the  deformation  is 
performed. 


2.2.  Extension  to  cone-beam  geometry 


Extension  of  the  motion  compensation  method  to  other  geometries  such  as  fan-beam  or  cone- 
beam  can  be  done  in  the  same  way  when  backprojection-type  algorithms  (Feldkamp  et  al  1984, 
Lange  and  Carson  1984)  are  used  for  reconstruction.  In  particular  for  the  circular  orbit  CBCT, 
the  Feldkamp  algorithm  commonly  built  in  commercial  machines  is  modified  to  accommodate 
the  motion  effects. 

Following  the  notation  in  Wang  et  al  (  1993),  with  R(6,  p,  g)  denoting  the  cone-beam 
projection  of  an  object  g(x,  y,  z)  at  gantry  rotational  position  0  and  detector  bin  (p,  g),  the 
Feldkamp  algorithm  can  be  expressed  as  the  ‘weighted  filtered  backprojection’, 


g(x,y,z)  =  - 


R(0 ,  u ,  v)  = 


\L 

/: 


(p  -  s)2 
p 


R(0,  u,  v)  d 0 


-R(0,  p,  v)f(u  -  p)dp, 


(8) 


(9) 


i-oo  y-p2 + p + v2 

where  p  is  the  source  to  axis  distance  (SAD),  s  =  —  x  sin#  +  y  cos  6,  R (6,  u,  i;)  is  the  filtered 
projection,  and  ( u ,  v)  is  the  intersecting  point  on  the  detection  plane  of  the  ray  coming  from 
the  cone  vertex  through  the  reconstruction  point  (x,y,z),  and  /(.)  is  the  filter  function.  More 
details  can  be  found  in  the  generalized  Feldkamp  algorithm  of  Wang  et  al  (1993).  Again, 
equation  (8)  can  be  rewritten  in  a  similar  form  to  (6),  and  we  can  apply  the  deformation  to  the 
single-angle  backprojected  objects  before  the  summation  operation  and  obtain  the  following 
approximate  formula: 

„2 


ge,(x,y,z)  = 
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R(9,,  u,  v), 


(10) 
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2.3.  Motion  model  derivation 

In  radiotherapy  clinics,  the  CBCT  on-board  imager  aims  to  accurately  position  the  patient 
in  the  treatment  room.  For  a  tumour  under  the  influence  of  respiration  motion,  4D  CT  has 
been  adopted  to  extract  organ  motion  information  for  4D  treatment  planning  (Keall  et  al  2005, 
Trofimov  et  al  2005,  Webb  2005).  The  4D  patient  model  derived  from  the  4D  CT  data  can  be 
applied  to  correct  motion  artefacts  in  CBCT  imaging. 

A  deformable  image  registration  model  is  required  to  obtain  the  4D  patient  model  from 
the  4D  CT  data.  For  convenience,  we  adopted  the  free-form  Spline  model  (B  Spline)  (Mattes 
et  al  2003)  in  this  work  to  register  the  4D  CT  images.  The  simplicity  and  yet  accuracy  of 
the  B Spline  method  make  it  a  preferred  tool  for  many  clinical  applications  (Lian  et  al  2004, 
Rohlfing  et  al  2003,  Rueckert  et  al  1999).  In  this  model,  a  lattice  of  user-defined  nodes 
is  overlaid  on  the  image.  Each  node  contains  a  deformation  vector,  whose  components  are 
determined  by  an  optimization  procedure.  The  deformation  at  any  point  of  the  image  is 
calculated  by  Spline  interpolation  of  the  adjacent  node  values.  One  advantage  of  the  B Spline 
model  is  that  the  nodes  are  locally  controlled,  such  that  the  displacement  of  an  interpolation 
point  is  influenced  only  by  the  nearest  grid  points  and  changing  a  lattice  node  only  affects  the 
transformation  regionally,  making  it  efficient  in  describing  local  deformations.  Suitable  node 
deformations  are  obtained  using  the  gradient-based  L-BFGS  algorithm  (Liu  and  Nocedal  1989, 
Schreibmann  and  Xing  2005),  which  iteratively  varies  the  deformations  until  the  registration 
metric,  a  mathematical  measure  of  similarity  between  images,  is  minimized.  The  normalized 
cross  correlation  metric  was  used  as  the  metric  in  this  study.  With  deformable  registration,  all 
phases  can  be  registered  to  one  particular  phase,  for  example  at  t  =  0,  resulting  in  a  series  of 
transformations  representing  a  temporal  sequence  of  3D  deformation  fields,  or  in  other  words, 
a  4D  model  of  organ  motion. 

2.4.  Computer  simulations 

The  performance  of  the  proposed  algorithm  was  tested  by  computer  simulations  in  this  work 
for  both  in-plane  motion  and  3D  motion.  First,  the  in-plane  motion  was  studied  with  a 
2D  digital  dynamic  phantom,  capable  of  both  rigid  translational  motion  and  non-rigid  shape 
change.  As  shown  in  figure  1,  the  phantom  consisted  of  several  ellipses  of  different  sizes  and 
densities.  The  whole  phantom  varied  its  shape  constantly  during  the  data  acquisition,  and  an 
internal  circle  was  also  moving  vertically  with  an  amplitude  of  1 .5  cm.  The  motion  period  was 
4.0  s,  of  which  the  inhale  stage  took  3.0  s  and  the  exhale  stage  took  1.0  s.  The  phantom  size 
was  512  x  512  pixels,  with  the  pixel  size  of  1  mm.  The  projections  were  simulated  by  line 
integrals  of  the  phantom  density  with  parallel-beam  geometry.  For  the  slow-rotating  CT,  the 
simulated  360°  scan  took  40.0  s  and  generated  512  equally  time-spaced  projections.  In  other 
words,  the  projections  were  acquired  at  time  t  =  i  x  40/511  s,  /  =0,1, ...,511,  with  the 
phantom  moving  correspondingly.  The  time  it  took  for  acquiring  each  single  projection  was 
neglected.  For  comparison,  a  simulated  full  scan  of  the  stationary  phantom  characterizing  its 
state  at  t  =  0  was  also  performed. 

Furthermore,  Gaussian  noise  was  added  into  the  imaging  process  to  test  the  robustness 
of  the  proposed  motion  correction  strategy.  Specifically,  30%  Gaussian  noise  was  added 
into  each  phase  of  the  motion  phantom  as  well  as  the  simulated  projections.  Therefore,  the 
deformation  fields  (or  motion  model)  derived  from  the  phantoms  were  subjected  to  the  impact 
of  the  noise,  and  the  errors  of  the  derived  motion  model  will  subsequently  propagate  into  the 
image  reconstruction  process. 

Finally,  3D  circular  orbit  CBCT  was  simulated  with  a  more  realistic  digital  phantom  to  test 
the  motion  correction  method.  In  this  study,  the  phantom  was  constructed  with  4D  CT  images 
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Figure  1.  A  2D  dynamic  phantom  for  CT  simulations  with  a  complex  mode.  The  whole  phantom 
changes  its  shape  as  indicated  by  the  dotted  line  with  a  period  of  3.52  s.  The  volume  of  the  phantom 
was  kept  constant  during  the  movement.  The  inner  right  circle  moves  vertically  with  the  same 
period  (plot  was  not  drawn  to  a  real  scale). 


Figure  2.  A  3D  deformable  phantom  for  CBCT  simulations.  The  phantom  was  constructed  from 
4D  CT  images  of  a  lung  patient.  The  CT  numbers  were  converted  to  the  attenuation  coefficient 
and  the  top  row  is  the  volume  rendering  of  the  data.  From  left  to  right  it  contains  the  transaxial, 
coronal  and  sagittal  images  of  the  phantom. 


of  a  lung  patient.  The  CT  numbers  (Hounsfield  unit)  were  converted  into  the  attenuation 
coefficient  and  a  set  of  eight  volumetric  data  corresponding  to  eight  respiratory  phases  of  the 
patient  were  obtained,  each  with  an  array  size  of  256  x  256  x  64  voxels.  Figure  2  shows  the 
phantom  at  the  first  phase  with  three  views  and  with  volume  rendering.  In  the  simulations, 
the  cone-beam  source  to  detector  distance  was  1500  mm,  and  1000  mm  to  the  centre  of  rotation. 
The  detector  size  was  512  mm  x  128  mm.  A  total  of  256  projections  were  simulated,  with  the 
phantom  moving  from  one  phase  to  another  at  subsequent  projection  angles.  For  comparison, 
a  static  CBCT  simulation  was  also  performed  with  the  same  phantom  for  the  first  phase.  The 
free-form  B  Spline  deformable  registration  was  applied  to  the  eight  phases  of  the  phantom  to 
extract  the  motion  model,  and  the  projection  data  were  then  reconstructed  with  the  conventional 
and  modified  Feldkamp  algorithm. 
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Figure  3.  Comparison  of  the  sinogram  with  moving  phantom  and  stationary  phantom.  Left: 
simulated  projections  with  the  2D  dynamic  phantom;  right:  simulated  projections  with  the  same 
phantom  with  the  motion  ‘switched  off’. 


3.  Results 

3.1.  2D  simulations 

Figure  3  shows  the  simulated  motion  corrupted  projections  with  the  2D  dynamic  phantom,  as 
well  as  the  projections  of  the  same  phantom  with  motion  ‘switched  off’.  The  inconsistency 
in  the  sinogram  of  the  moving  phantom  is  clearly  seen  in  the  picture.  To  compensate  for 
the  motion  in  the  image  reconstruction,  we  first  derived  the  motion  model  for  the  dynamic 
phantom.  The  deformation  field  at  each  phase  with  respect  to  the  first  phase  was  obtained  by 
the  B Spline  deformable  registration  described  earlier.  As  an  example,  we  illustrate  in  figure  4 
the  deformation  fields  of  phase  i  =  10  relative  to  phase  i  =  0  at  each  region  of  the  phantom. 
The  arrows  in  the  figure  show  the  direction  of  the  movement  at  each  pixel  and  the  colour 
indicates  the  amplitude  of  the  deformation  vectors. 

By  applying  the  deformation  in  the  backprojection  step  of  the  reconstruction,  the  images 
with  reduced  motion  artefacts  can  be  obtained.  The  reconstructed  images  with  and  without 
motion  correction  are  shown  in  figure  5,  in  which  the  left  one  is  the  phantom,  the  middle 
is  the  reconstructed  image  without  motion  correction,  and  the  right  one  is  the  reconstructed 
image  with  motion  correction.  It  is  observed  that  the  outer  boundary  of  the  motion  phantom  is 
corrected  with  the  approach,  and  the  shape  of  the  small  moving  circle  is  also  restored.  Figure  6 
shows  the  vertical  profiles  through  the  moving  circle,  in  which  it  is  found  that  the  intensity  of 
the  moving  circle  is  in  accordance  with  the  phantom  value. 

In  order  to  test  the  robustness  of  the  proposed  method,  we  further  added  30%  Gaussian 
noise  into  the  simulation  to  create  uncertainties  in  both  projections  and  the  derived  deformation 
fields.  Although  the  deformable  registration  was  affected  by  the  noise  (see  figures  7  and  4),  it 
is  found  that  the  density  of  the  moving  part  of  the  phantom  was  recovered  correctly  and  the 
outer  boundary  of  the  phantom  was  corrected  as  well,  see  figure  8. 

3.2.  3D  cone-beam  simulations 

The  proposed  motion  compensation  method  was  further  tested  under  cone-beam  geometry  with 
a  4D  anthropomorphic  digital  phantom.  Figure  9  shows  the  eight  phases  of  the  phantom.  The 
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Figure  4.  Deformation  field  obtained  by  B  Spline  registration  of  phase  i  =  10  and  phase  i  =  0  for 
the  motion  phantom  shown  in  figure  1 . 


Figure  5.  Phantom  and  reconstructed  images.  Left:  phantom  image;  middle:  reconstructed  image 
without  correction;  right:  reconstructed  image  with  motion  correction. 


derived  deformation  field  for  phase  1  to  phase  0  is  illustrated  in  figure  10.  The  reconstructed 
images  with  the  conventional  and  modified  Feldkamp  algorithm  are  presented  in  figure  1 1 , 
in  which  the  first  column  shows  the  reconstructed  images  of  CBCT  projection  of  the  moving 
phantom  with  the  conventional  Feldkamp  algorithm,  and  the  middle  column  shows  the  same 
data  reconstructed  with  the  proposed  motion  compensation  algorithm.  For  comparison, 
the  static  simulation  with  the  phantom  at  phase  0  was  reconstructed  with  the  conventional 
Feldkamp  algorithm  and  is  shown  in  the  last  column.  From  top  to  bottom  in  the  figure  are 
the  axial,  coronal  and  sagittal  views,  and  the  last  row  shows  zoom-in  images  of  the  region  of 
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Figure  6.  Vertical  profiles  through  the  inner  moving  circle  for  images  shown  in  figure  5. 


Figure  7.  Deformation  field  obtained  in  the  presence  of  Gaussian  noise  by  B Spline  registration 
for  phase  i  =  10  and  phase  i  =  0  for  the  2D  dynamic  phantom. 


interest  illustrated  in  the  axial  view.  It  is  found  that  the  motion  artefacts  were  reduced  and  the 
tumour  intensity  and  shape  were  restored  by  our  proposed  method. 
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Figure  8.  Reconstructed  images  with  noise  present.  Left:  reconstructed  image  without  correction; 
right:  reconstructed  image  with  motion  correction.  The  derived  motion  model  was  affected  by  the 
noise;  however,  the  reduction  of  the  tumour  distortion  was  still  observed. 


4.  Discussion  and  conclusion 

There  are  generally  two  approaches  to  the  removal  of  organ  motion  artefacts  in  CT  imaging. 
The  first  is  to  sort  the  projections  according  to  their  phases  and  then  reconstruct  the  sorted 
projections  separately  to  obtain  time-resolved  CT  images.  Another  approach  is  to  incorporate 
the  patient’s  motion  model,  which  describes  how  each  point  of  the  object  moves  during 
the  data  acquisition,  into  the  image  reconstruction  process.  For  CBCT,  the  first  approach 
is  less  adequate  because  there  are  not  sufficient  projects  for  a  given  phase  to  warrant  high 
quality  4D  CBCT  images,  unless  multiple  gantry  rotations  are  performed.  In  this  work,  we 
have  demonstrated  the  feasibility  of  incorporating  a  4D  motion  model  into  the  CBCT  image 
reconstruction.  The  motion  model  used  in  this  work  is  derived  on  a  patient- specific  basis, 
therefore  could  be  more  suitable  than  other  mathematical  models  for  a  patient  study.  In 
radiotherapy  clinics,  4D  CT  has  been  adopted  to  identify  the  tumour  motion  and  to  guide  the 
4D  treatment  planning  process.  Using  the  4D  CT  data  for  the  removal  of  motion  artefacts  in 
CBCT  acquisition  seems  to  be  logical  and  fits  well  with  the  data  flow  of  modern  image-guide 
radiotherapy  (IGRT). 

Incorporating  a  motion  model  into  the  de-convolution  process  is  a  general  approach  in 
estimation  theory  and  has  been  applied  to  CT  and  MRI  reconstruction  and  radiation  dose 
optimization  (Li  and  Xing  2000,  Pugachev  and  Xing  2002).  In  this  process,  the  motion  model 
serves  as  a  priori  knowledge  of  the  system  and  provides  valuable  partial  guidance  to  the 
searching  process.  It  should  be  noted  that  a  potential  shortcoming  of  CBCT  reconstruction 
with  inclusion  of  the  4D  patient  model  derived  at  an  earlier  time  point  is  that  it  does  not 
consider  possible  change  in  the  patient’s  organ  motion  pattern.  Indeed,  the  motion  pattern 
may  change  in  both  spatial  and  temporal  domains  to  some  degree  from  the  planning  CT  when 
the  treatment  CBCT  is  acquired.  This  represents  an  important  challenge  not  only  to  the  4D 
CBCT  reconstruction  proposed  here,  but  also  to  the  implementation  of  4D  radiation  therapy 
(e.g.,  tumour  tracking  based  on  4D  planning).  Other  uncertainties  of  the  motion  model  include 
spatial  and  temporal  errors,  for  example,  due  to  accuracy  in  aligning  the  CBCT  coordinate 
system  with  the  4D  CT  coordinate  system  and  in  determining  the  phase  of  each  individual 
CBCT  projection  or  4D  CT  images.  However,  unless  the  organ  motion  pattern  is  completely 
reversed,  we  anticipate  that  the  proposed  technique  will  help  to  reduce  the  motion  artefacts 
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Figure  9.  The  eight  phases  of  the  4D  anthropomorphic  phantom. 


even  if  a  slight  change  in  the  motion  pattern  occurs.  In  practice,  a  consistent  respiratory 
pattern  should  be  maintained  through,  for  example,  a  proper  voice  coaching,  during  the  CBCT 
data  acquisition  process.  This  research  also  highlights  the  importance  of  the  development  of 
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Figure  10.  The  derived  deformation  fields  of  the  3D  phantom  at  phase  1  relative  to  phase  0  by 
B  Spline  deformable  registration. 


a  robust  organ  motion  tracking  system  in  the  future.  With  this,  one  may  obtain  the  4D  patient 
model  during  the  4D  CBCT  process  and  then  incorporate  it  into  the  CBCT  image  reconstruction 
process.  A  study  to  systematically  examine  the  sensitivity  of  the  4D  CBCT  reconstruction 
against  various  possible  variations  of  the  patient  motion  model  should  be  of  scientific  and 
practical  interest.  Because  of  the  large  scope  of  the  work,  we  focus  this  manuscript  on  setting 
up  the  theoretical  framework  of  4D  CBCT  reconstruction  and  leave  the  sensitivity  study  for 
the  future  investigations. 

In  the  proposed  approach,  when  deformation  field  is  applied  with  backprojection,  the 
reconstructed  image  is  locally  corrected,  which  is  same  as  Rietch  et  al  did  with  the  pixel- 
specific  ‘magnification  and  displacement’  motion  model.  It  should  be  noted  that  the  local 
correction  does  not  correct  the  artefacts  that  are  left  to  other  regions,  for  example,  regions  that 
have  no  motion  at  all  but  still  suffer  from  the  artefacts  produced  by  the  moving  part.  In  general, 
an  iterative  method  could  be  utilized  to  incorporate  the  motion  model  for  image  reconstruction 
with  better  motion  artefacts  correction.  The  ordered  subset  convex  (OSC)  algorithm  (Kole 
and  Beekman  2005,  Manglos  et  al  1992)  could  be  a  good  choice  for  the  CBCT  task,  since  the 
subset  can  be  naturally  correlated  to  the  phase  set  for  an  efficient  implementation.  However, 
the  presented  approach  should  be  much  faster  than  iterative  methods  with  sufficient  restoration 
of  the  moving  target  information. 

In  summary,  we  have  proposed  a  motion  correction  method  for  reconstructing  CT  images 
under  the  influence  of  intra-fraction  organ  motion.  The  technique  is  conceptually  interesting 
and  may  find  a  natural  application  in  routine  CBCT-based  patient  positioning  for  improved 
accuracy.  The  computation  is  efficient,  involving  the  standard  filtered  backprojection  step 
and  deformation  for  each  projection  phase.  The  method  deals  with  not  only  in-plane  motion, 
but  also  3D  complex  motion  via  deformable  registration.  It  enhances  the  quality  of  the 
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Figure  11.  Reconstructed  images  for  the  CBCT  simulations.  Left  column:  reconstructed  images 
without  motion  correction;  middle  column:  images  with  motion  correction;  right  column:  images 
reconstructed  with  the  conventional  Feldkamp  algorithm  for  projections  with  stationary  phantom 
(the  3D  phantom  at  phase  0).  From  top  to  bottom  are  the  axial,  coronal  and  sagittal  views,  and  the 
last  row  contains  the  zoom-in  images  of  the  region  of  interest  shown  in  the  first  row. 


reconstructed  image,  corrects  for  the  density,  and  restores  the  shape  of  the  moving  tumour, 
therefore,  can  provide  better  target  localization. 
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Four-dimensional  (4D)  CT  is  useful  in  many  clinical  situations,  where  detailed  abdominal  and 
thoracic  imaging  is  needed  over  the  course  of  the  respiratory  cycle.  However,  it  usually  delivers  a 
larger  radiation  dose  than  the  standard  three-dimensional  (3D)  CT,  since  multiple  scans  at  each 
couch  position  are  required  in  order  to  provide  the  temporal  information.  Our  purpose  in  this  work 
is  to  develop  a  method  to  perform  4D  CT  scans  at  relatively  low  current,  hence  reducing  the 
radiation  exposure  of  the  patients.  To  deal  with  the  increased  statistical  noise  caused  by  the  low 
current,  we  proposed  a  novel  4D  penalized  weighted  least  square  (4D-PWLS)  smoothing  method, 
which  can  incorporate  both  spatial  and  phase  information.  The  4D  images  at  different  phases  were 
registered  to  the  same  phase  via  a  deformable  model,  thereby,  a  regularization  term  combining 
temporal  and  spatial  neighbors  can  be  designed  for  the  4D-PWLS  objective  function.  The  proposed 
method  was  tested  with  phantom  experiments  and  a  patient  study,  and  superior  noise  suppression 
and  resolution  preservation  were  observed.  A  quantitative  evaluation  of  the  benefit  of  the  proposed 
method  to  4D  radiotherapy  and  4D  PET/CT  imaging  are  under  investigation.  ©  2005  American 
Association  of  Physicists  in  Medicine.  [DOI:  10.1118/1.2122567] 

Key  words:  4D  CT,  low  dose,  4D  radiotherapy  planning,  PET/CT,  PWLS,  deformable  registration 


I.  INTRODUCTION 

In  radiation  therapy,  respiratory  motion  poses  significant 
challenges  for  tumors  in  the  thorax  or  abdomen.  It  can  distort 
the  shape  of  an  object,  degrade  the  anatomic  position  repro¬ 
ducibility  during  imaging,  and  necessitate  larger  margins 
during  radiotherapy  planning.  It  also  causes  inaccuracy  in 
estimating  the  tumor  volume,  thereby  preventing  an  effective 
dose  escalation  for  the  treatment  of  a  target  tumor.  These 
issues  make  it  difficult  to  achieve  the  desired  goals  of  con¬ 
formal  radiotherapy.  Four-dimensional  (4D)  CT  scans,  ac¬ 
quired  synchronously  with  a  respiratory  signal,  provide  not 
only  the  three-dimensional  (3D)  spatial  information,  but  also 
temporal  changes  of  the  anatomy  as  a  function  of  the  respi¬ 
ratory  phase  during  the  imaging,  and  can  therefore  be  em¬ 
ployed  in  4D  treatment  planning  to  explicitly  account  for  the 
respiratory  motion.1,2 

To  acquire  4D-CT  scans,  a  position-monitoring  system 
used  for  respiratory  motion  tracking  is  usually  interfaced 
with  the  CT  scanner,  so  that  the  CT  data  are  acquired  in 
correlation  with  real-time  positioning.  One  scanning  proto¬ 
col,  which  has  recently  been  developed  to  achieve  4D-CT 
imaging  with  the  multislice  CT  scanner,3,4  uses  a  “step-and- 
shoot”  technique  (cine  mode),  in  which  CT  projections  are 
acquired  over  a  complete  respiratory  cycle  at  each  couch 
position.  The  period  of  each  CT  acquisition  segment  is  time 
stamped  with  an  “x-ray  ON”  signal  and  recorded  by  the 
tracking  system.  The  4D-CT  data  are  subsequently  sorted 
into  groups  according  to  their  phases  of  the  breathing  cycle. 
Typical  parameters  for  thoracic  imaging  are  0.45  s  cine  in¬ 
tervals  for  a  duration  slightly  longer  (1  s)  than  a  full  respi¬ 
ratory  cycle  at  each  couch  position,  a  0.5  s  gantry  rotation, 


140  kVp,  175  mA,  a  2.5  mm  slice  thickness,  and  10  mm 
couch  increments  for  a  four-row  scanner.5  This  data  acquisi¬ 
tion  technique  takes,  on  average,  about  one  minute,  depend¬ 
ing  on  the  patient  breathing  period  and  axial  dimension  of 
the  scan  range. 

When  a  modern  multislice  CT  is  used  for  a  regular  clini¬ 
cal  exam,  the  dose  received  by  the  patient  may  approach 
10  mSv  for  the  head,  and  20  mSv  for  the  chest  or  abdomen. 
With  a  4D  acquisition,  since  a  patient  is  scanned  multiple 
times  at  each  couch  position  during  the  imaging,  the  radia¬ 
tion  exposure  will  be  considerably  higher  than  the  regular 
CT  scan  (up  to  one  order  of  magnitude  higher).  The  effective 
dose  reduction  is  thus  highly  desirable  for  a  clinical  applica¬ 
tion  of  the  cutting-edge  4D-CT  scanning  technology.  In  this 
work  we  develop  a  novel  method  to  perform  the  4D  scan  at 
a  lower  current  and  retrospectively  generate  images  compa¬ 
rable  to  a  high-mA  4D  scan,  hence  reducing  the  patient  ra¬ 
diation  dose.  The  central  idea  is  to  map  different  phases  into 
a  particular  phase  through  deformable  model  registration 
methods,  and  to  improve  the  image  for  this  phase  by  statis¬ 
tical  estimation  from  the  registered  images.  In  the  following, 
we  describe  the  method  and  the  results  of  phantom  and  pa¬ 
tient  studies  in  detail. 

II.  METHODS  AND  MATERIALS 
A.  Phantom  data  acquisition 

Two  phantom  experiments  were  carried  out  in  this  work: 
one  is  a  commercial  calibration  phantom  CatPhan®  600  (The 
Phantom  Laboratory,  Inc.,  Salem,  NY),  as  shown  in  Fig.  1, 
and  the  other  is  an  anthropomorphic  thorax  phantom  (Fig.  2). 
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CTP404,  HU  calibration,  spatial  linearity 
CTP591,  Bead  geometry 
CTP528,  21  line  pair  high  resolution 
CTP528.  Point  source 

CTP515,  Subslice  and  supra-slice  low  contrast 
CTP486,  Solid  image  uniformity  module 


Distance  from  section  1  center 
Omm 
32.5mm 
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Fig.  1.  CatPhan®  600  phantom  and  its  module  specification.  The  module 
section  CTP528  has  high  resolution  multiple  line  pairs  and  was  used  in  this 
work  to  test  the  spatial  resolution  loss  during  the  proposed  image 
processing. 


To  acquire  the  4D  CT  data,  each  phantom  was  placed  on  the 
top  of  a  platform  capable  of  sinusoidal  motion  along  the 
cranial-caudal  direction.  The  amplitude  of  motion  varied  dis¬ 
cretely  from  1  to  6  cm,  and  the  period  was  continuously  ad¬ 
justable  from  0.5  s  to  1  min.  A  Varian  Real-time  Position 
Management  (RPM)  respiratory  gating  system  (Varian  Medi¬ 
cal  Systems,  Palo  Alto,  CA)  was  used  to  record  the  motion 
by  tracking  two  infrared  reflective  markers,  rigidly  mounted 
on  a  plastic  block  on  the  top  of  the  phantom,  by  means  of  an 
infrared  video  camera  mounted  on  the  PET/CT  table  (see 
Fig.  3).  The  motion  amplitude  was  displayed  as  a  function  of 
time  at  a  rate  of  30  Hz  on  the  RPM  workstation,  and  the  data 
were  recorded  for  the  entire  4D-CT  scan.  During  the  scan, 
the  RPM  system  also  recorded  the  state  of  a  TTL  (Transistor- 
Transistor  Logic)  “x-ray  ON”  signal  from  the  CT  scanner, 
indicating  the  acquisition  time  of  the  CT  data,  thereby  time 
stamping  each  CT  slice  with  respiratory  motion.  Figure  4 
shows  one  example  of  the  motion  wave  recorded  by  the 


Fig.  2.  The  experimental  design  showing  the  power  supply  that  drives  the 
electric  motor,  which  in  turn  moves  the  motion  platform  and  phantom  sinu¬ 
soidally  in  the  cranial-caudal  direction  as  it  moves  through  the  CT  bore. 


Fig.  3.  The  RPM  infrared  camera  and  illuminator  system  is  mounted  at  the 
foot  of  the  couch  and  used  to  track  the  motion  of  infrared  reflectors  placed 
on  the  top  of  the  thorax  phantom  or  the  diaphragm  of  the  patient. 


RPM  system,  in  which  the  diamond  markers  indicate  the 
motion  status  when  a  CT  set  is  generated  at  each  couch  po¬ 
sition. 

A  Discovery  ST  PET/CT  scanner  (Discovery  ST/ 
LightSpeed  8-slice;  General  Electric  Medical  Systems)  was 
used  in  this  study.  For  CT  scans,  it  has  a  50  cm  transaxial 
field  of  view  and  an  available  slice  thickness  of  between 
0.625  and  10.0  mm.  The  tube  current  can  be  varied  between 
10  and  440  mA,  and  the  available  tube  voltage  settings  are 
80,  100,  120,  and  140  kVp.  The  detector  coverage  in  the 
cranial-caudal  direction  is  2  cm,  and  the  minimum  gantry 
rotation  time  is  0.5  s. 

In  our  studies,  both  phantoms  moved  with  a  period  of 
5.0  s  and  amplitude  of  2  cm.  The  CT  data  were  acquired  in 
cine  mode,  at  120  kVp,  10  mA,  with  cine  CT  axial  field  of 
view  of  20  mm  (8  X  2.5  mm  slice  thickness).  The  x-ray  tube 
angular  velocity  was  set  to  0.5  s/rotation.  The  “shoot”  period 
at  each  couch  position  (cine  duration)  was  set  slightly  longer 
than  the  motion  period  to  6.1  s,  and  the  cine  interval  between 
images  was  0.45  s.  Each  image  was  reconstructed  with  360° 
of  data.  The  scan  covered  10  cm  (5  couch  positions)  with  the 
total  acquisition  time  of  about  30  s.  For  CatPhan®  600,  the 
image  can  be  compared  with  the  specification  of  the  phantom 
manufacture;  to  compare  the  image  quality  for  the  thorax 
phantom,  another  4D  scan  was  repeated  at  100  mA  with  all 
other  parameters  kept  the  same. 

B.  Patient  data  acquisition 

A  clinical  4D-CT  patient  study  was  performed  on  the 
same  scanner.  The  patient  was  asked  to  breathe  normally 
during  the  scan.  The  plastic  block  with  two  infrared  reflec¬ 
tive  markers  was  taped  on  the  top  of  the  patient’s  abdomen, 
placed  medially  and  a  few  cm  inferior  to  the  xiphoid  pro¬ 
cesses.  The  respiratory  signal  of  the  patient  from  the  RPM 
was  recorded  and  synchronized  with  CT  data  acquisition, 
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Fig.  4.  An  illustration  of  the  respiratory  motion  wave¬ 
form  recorded  by  the  RPM.  Each  diamond  marker  re¬ 
flects  the  motion  status  when  a  CT  volumetric  image 
(for  one  couch  position)  is  generated.  The  example  here 
shows  that  16  CT  images  were  generated  in  6.1  s  at 
each  couch  position. 


similar  to  the  phantom  scan.  The  axial  coverage  of  the  pa¬ 
tient  was  25  cm.  Other  scan  parameters  were  90  mA, 
120  kV,  and  2.5  mm  slice  thickness.  The  dose  to  the  patient 
from  the  4D-CT  was  80.2  mGy.  The  cine  interval  between 
images  was  at  0.45  s.  Each  image  reconstruction  used  360° 
of  data  corresponding  to  0.8  s  duration. 

C.  Deformable  registration  for  phase  images 

After  the  scan  data  were  prospectively  reconstructed  at 
the  PET/CT  scanner,  both  the  CT  images  and  the  correspond¬ 
ing  motion  data  recorded  by  the  RPM  system  were  trans¬ 
ferred  to  a  GE  Advantage  Workstation  (GE  Medical  Systems, 
Waukesha,  WI).  The  “Advantage  4D”  software  on  the  work¬ 
station  simultaneously  displays  the  CT  images  and  the  mo¬ 
tion  data,  and  sorts  the  cine  images  into  a  set  of  respiratory 
phase  images.  For  both  the  phantom  and  the  patient  scans  in 
this  study,  a  total  of  10  phases  were  created,  with  phase 
intervals  of  10%  of  the  respiratory  cycle. 

The  phase  images  obtained  from  low-mA  scans  are  quite 
noisy  (see  examples  in  Figs.  6  and  8,  later)  and  usually  can¬ 
not  be  used  directly  for  4D  radiation  treatment  planning.  Be¬ 
fore  applying  our  4D  image  enhancement  technique  to  these 
noisy  phase  images,  a  registration  step  was  performed  in 
order  for  the  temporal  information  to  be  appropriately  incor¬ 
porated  in  the  next  step.  To  account  for  the  complex  organ 
motion,  we  have  previously  investigated  a  few  deformable 
models,6,7  and  adopted  the  free-form  spline  model  (B Spline) 
in  this  work.8  Its  simplicity  and  yet  accuracy  make  it  a  pre¬ 
ferred  tool  for  clinical  applications.9-11  In  this  model,  a  lat¬ 
tice  of  user-defined  nodes  is  overlaid  on  the  image.  Each 
node  contains  a  deformation  vector,  whose  components  are 
to  be  determined  by  the  optimization  procedure.  The  defor¬ 
mation  at  any  point  of  the  image  is  calculated  by  spline 
interpolation  of  closest  nodes  values.  Unlike  other  spline 
models,  the  B Splines  are  locally  controlled,  such  that  the 
displacement  of  an  interpolation  point  is  influenced  only  by 
the  closest  grid  points  and  changing  a  lattice  node  only  af¬ 
fects  the  transformation  regionally,  making  it  efficient  in  de¬ 
scribing  local  deformations.  Suitable  node  deformations  are 
solved  using  the  gradient-based  algorithm  L-BFGS7  12  due  to 
its  superior  performance  in  large-scale  optimization  prob¬ 
lems.  The  optimizer  iteratively  varies  deformation  values  to 
minimize  the  metric,  a  mathematical  measure  of  similarity 
between  images.  The  normalized  cross  correlation  (NCC) 
metric  (described  in  the  next  section)  was  used  here  since  the 
registration  was  applied  on  images  acquired  under  identical 
settings.13 


To  improve  the  4D  images,  for  example,  at  phase  i,  all 
other  phases  were  registered  to  phase  i  and  resulted  in  trans¬ 
formations  representing  a  temporal  sequence  of  3D  deforma¬ 
tion  fields,  or  in  other  words,  a  4D  model  of  organ  motion. 
Each  phase  throughout  the  respiratory  cycle  can  be  subse¬ 
quently  warped  with  the  deformation  field  to  map  the  points 
in  each  phase  to  the  corresponding  points  in  the  reference 
image,  i.e.,  the  phase  i  image. 

D.  4D-CT  image  enhancement 

The  goal  of  image  enhancement  is  essentially  to  obtain  a 
more  accurate  estimation  of  the  “true”  intensity  (CT  number) 
of  each  voxel  in  the  image  being  processed.  In  our  case,  the 
information  available  here  included  the  image  at  phase  i,  and 
nine  other  images  that  were  previously  registered  to  phase  i. 
Ideally,  if  the  nine  phases  are  perfectly  registered  to  phase  /, 
then  each  of  them  can  be  considered  a  repeated  scan  for 
phase  i.  Therefore,  it  is  natural  to  average  them  with  phase  i 
to  get  the  improved  image.  In  reality,  however,  there  are 
always  errors  in  registration,  and  not  all  voxels  will  be  able 
to  match  exactly.  Simple  averaging  will  thus  lead  to  blurring 
in  the  image.  In  the  following,  we  describe  a  statistical 
method  to  incorporate  the  available  information  into  a  penal¬ 
ized  weighted  least  square  (PWLS)  objective  function  to 
achieve  the  optimal  estimation  of  the  true  image  at  phase  i. 
The  method  extends  the  conventional  PWLS  method14,15  into 
four  dimensions  to  include  the  temporal  information,  and 
from  here  on  we  refer  to  it  as  the  4D-PWLS  smoothing 
method. 

Let  X  denote  the  “4D  deformed  volumetric  image”  X 

-  (-U  •••  ^N,XN+1,XN+2,  ...  ,X2n,  •••  >*(M-l)V+l>*(M-l)V+2> 

...  ,xMN)',  where  M  is  the  total  number  of  phases  after  de¬ 
formable  registrations,  N  is  the  total  number  of  voxels  in 
each  phase  image,  and  the  prime  4  denotes  the  transpose 
operation.  Let  /i  be  the  unknown  true  4D  image 

A*  =  ,  yL62  ,  .  .  .  ,  /%,  /%+1 ,  /%+2,  •  •  •  >  M2AG  •  •  •  >  M(M-1)V+1  > 

/^(M-\)N+2y  •  •  •  ■  The  4D-PWLS  solution  is  to  find  jd 

that  minimizes  the  following  objective  function: 

$(/*)  =  I(X  -  At)'2-‘(X  -  At)  +  aH(fi) 

At  =  argmin<3>(At),  (1) 

where  the  first  term  is  the  weighted-least- square  part,  in 
which  the  variance-covariance  matrix  2  is  diagonal  with  the 
jth  entry  oj,  j=  1 ,2, . . . ,  AX  M,  if  the  intensity  measurement 
at  each  voxel  of  the  4D  image  is  statistically  independent.  To 
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(a)  (b) 

Fig.  5.  An  illustration  of  the  registration  error,  (a)  A  perfect  match  of  the 
corresponding  voxels;  (b)  when  errors  exist,  the  voxels  transformed  by  reg¬ 
istration  may  not  match  exactly  their  targeted  locations. 


provide  a  relatively  accurate  estimation  of  the  variance  oj  at 
each  voxel  j,  a  practical  way  is  to  select  a  “window”  for  each 
voxel  and  calculate  the  sample  variance  in  the  window  as  the 
estimation.  For  4D  CT,  the  window  can  be  selected  as  a  set 
that  includes  neighbors  in  both  the  spatial  and  temporal  do¬ 
mains,  for  example,  the  four  nearest  spatial  neighbors  within 
frame  i,  and  (M- 1)  neighbors  at  the  corresponding  locations 
in  the  registered  frames. 

The  second  term  in  (1)  is  the  regularizing  penalty  term 
that  encourages  the  neighborhood  smoothness,  and  the  coef¬ 
ficient  a>0  is  a  constant  to  weight  the  penalty.16,17  Because 
both  spatial  and  temporal  neighbors  are  involved,  in  this 
work,  a  simple  quadratic  form  was  selected  to  penalize  the 
two  kinds  of  neighbors: 


H(a*)  =  ^2  2  WjXinj- nk)2, 


9  ““  "  JKrs  ' 

Z  j  k^Nj  Z 


(2) 


where  Nj  is  the  set  of  “4D  neighbors”  of  the  jth  pixel,  which 
include  eight  in-plane  spatial  neighbors  and  nine  temporal 
(phase)  neighbors  in  our  studies.  As  illustrated  in  Fig.  5,  if 
each  phase  is  perfectly  registered  to  a  reference  image  (for 
example,  phase  /),  then  the  corresponding  voxels  overlap 
with  each  other  (Fig.  5(a));  when  errors  are  present  in  regis¬ 
tration  (due  to  noise,  deformable  model,  interpolation,  etc.), 
the  points  in  the  registered  image  may  not  be  transformed  to 
the  corresponding  point  in  the  reference  image,  but  some 
place  nearby,  as  shown  in  Fig.  5(b).  Therefore,  these  corre¬ 
sponding  voxels  in  the  registered  phases  can  be  considered 
as  the  new  neighbors  of  the  reference  image  of  phase  i,  in 
addition  to  its  own  3D  spatial  neighbors.  The  weights  wjk 
equal  1  for  horizontal  and  vertical  spatial  neighbors,  and 
1/ V2  for  diagonal  neighbors.  For  the  temporal  neighbors,  the 
weights  were  proportional  to  the  NCC  between  current  frame 
m  and  frame  i,  which  is  defined  as 


fir 


Jj=  1 


XM  m-  \)Nxj+(i-\)N 


2  =1  (*j+(m-l)Ar)22  j=1  (■ Xj+(i-l)N )' 


(3) 


The  proportion  is  determined  by  the  ratio  of  the  NCCs  in 
temporal  and  spatial  directions.  Specifically,  the  NCCs  of 


Fig.  6.  Phantom  study  for  the  4D-PWLS  method  with  the  CatPhan®  600 
phantom.  The  CTP528  section  (high-resolution  line  pairs)  is  compared  in 
two  window  width/level  settings  for  10  mA  CT  image  before  (top  row)  and 
after  (bottom  row)  the  proposed  4D-PWLS  smoothing.  The  left  column  is 
displayed  with  window  width  200,  level  80;  the  right  column  is  displayed 
with  window  width  1500,  level  150. 


adjacent  rows  in  frame  i  is  calculated  similarly  as  (3),  and  an 
averaged  similarity  metric  fs  is  obtained  for  frame  i  for  the 
spatial  direction.  Naively,  larger  the  similarity  measure  is 
(maximum  1),  higher  the  weight  should  be.  Therefore,  we  set 
the  weights  of  temporal  neighbors  to  be  fim/fs.  Note  that  for 
the  objective  function  <F(/n)  to  be  convex,  the  temporal 
weights  must  also  be  positive.  A  unique  simple  closed-form 
solution  to  (1)  exists: 

yw  =  (2_1  +  afl)-12-1X.  (4) 

Since  the  matrix  H  is  multidiagonal,  the  inversion  involved 
in  (4)  is  not  straightforward.  In  this  work  we  adopt  a  simple 
iterative  algorithm,  known  as  the  iterated  conditional  mode 
(ICM),18 


which  converges  to  the  unique  solution. 

III.  RESULTS 
A.  Phantom  study 

The  CatPhan®  600  phantom  is  capable  of  measuring  the 
spatial  resolution.  In  Fig.  6  we  compared  the  images  for  one 
slice  of  the  CatPhan®  600  that  contains  multiple  strips  on  a 
circle,  of  which  the  resolution  sections  are  ranging  from  1  to 
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Table  I.  Mean  and  standard  deviation  (SD)  of  the  CT  number  in  the  se¬ 
lected  ROIs  before  and  after  the  proposed  4D  smoothing  for  the  10  mA  scan 
of  the  CatPhan  600  phantom.  ROI-1  is  the  10  X  10  pixel  square  region  in  the 
center,  ROI-2  is  the  10  X  10  pixel  square  region  5.0  cm  above  the  center, 
ROI-3  is  the  10  X  10  pixel  square  region  8.5  cm  above  the  center. 


ROI-1 

ROI-2 

ROI-3 

Mean  CT  number  (HU) 
in  original  image 

94.75 

100.19 

12.5 

Mean  CT  number  after 

4D  smoothing  (HU) 

95.94 

100.12 

11.43 

SD  CT  number  (HU) 
in  original  image 

17.39 

15.76 

13.12 

SD  CT  number  after 

4D  smoothing  (HU) 

7.46 

6.24 

5.31 

21  lines  pairs  per  cm.  The  original  and  smoothed  CT  images 
are  compared  in  two  different  window  width/level  settings, 
in  which  the  left  column  is  displayed  at  window  width  200, 
level  80,  and  the  right  column  is  displayed  at  window  width 
1500,  level  150.  The  top  row  is  the  original  4D  CT  images 
acquired  at  10  mA,  and  the  bottom  is  the  smoothed  images 
with  a- 1.0  X  10-6.  Since  the  proposed  weighted  least  square 
algorithm  is  adaptive  to  the  local  variance  at  each  pixel  of 
the  image,  the  resolution  and  noise  in  the  smoothed  image 
could  be  variable  at  different  locations.  We  examined  three 
regions  of  interest  (ROIs)  for  measurements  of  the  average 
CT  number  and  standard  deviation,  each  containing  10 
X  10  pixels  in  a  uniform  rectangular  region.  The  first  ROI 
was  located  at  the  center  of  the  image,  ROI-2  was  5.0  cm 
above  the  center,  and  ROI-3  was  8.5  cm  above  the  center. 
The  results  are  listed  in  Table  I.  It  is  found  that  the  smoothed 
image  had  the  similar  average  CT  number  as  the  original 
image  for  the  measured  ROIs,  while  the  standard  deviation 
of  the  CT  numbers  were  reduced  from  17.39,  15.76,  13.12  to 
7.46,  6.24,  and  5.31,  respectively.  Thus,  about  a  2.5X  reduc¬ 
tion  of  the  variation  of  CT  number  was  achieved.  Obviously, 
the  reduction  of  noise  itself  is  not  enough  to  assess  the  over¬ 
all  image  quality.  We  therefore  further  compared  the  two  sets 


(HU) 


Fig.  7.  Profiles  across  the  strips  indicated  by  the  squares  in  the  images  of 
Fig.  6.  The  dotted  line  shows  the  profile  of  the  original  10  mA  data,  and 
solid  line  is  the  profile  for  the  smoothed  data. 


of  images  for  the  spatial  resolution.  The  vertical  profiles 
across  the  strips  in  the  middle  right  (see  Fig.  6)  were  plotted 
in  Fig.  7,  where  only  a  small  resolution  loss  was  observed. 
The  four  peak  values  (or  the  “signals”)  from  left  to  right 
dropped  less  than  6%  from  1155,  1084,  1129,  1105  to 
1086.9,  1029.2,  1079.5,  1069.8,  respectively,  after  process¬ 
ing. 

In  the  human  thorax  phantom  study,  we  compared  the 
smoothed  CT  images  with  the  original  low-mA  (10  mA)  CT 
images,  as  well  as  the  original  high-mA  (100  mA)  CT  im¬ 
ages.  Figure  8  shows  their  coronal  images  at  five  different 
phases  of  0%,  20%,  40%,  60%,  and  80%.  The  left  column 
contains  the  4D  CT  images  acquired  at  100  mA  after  sorting 
by  the  GE  4D  Advantage  software.  The  middle  column 
shows  the  images  at  the  same  phases  acquired  at  10  mA,  and 
the  4D-PWLS  smoothed  images  with  a- 2. OX  10-6  based  on 
the  low-current  CT  data  are  displayed  in  the  right  column. 
The  effective  reduction  of  the  image  noise  can  be  observed  at 
a  certain  cost  of  the  spatial  resolution.  To  quantify  this,  we 


Fig.  8.  Phantom  study  for  the  4D-PWLS  method  with 
the  thorax  phantom.  The  left  and  middle  columns  are 
the  original  phases  obtained  from  the  GE  Advantage 
Workstation,  for  100  and  10  mA,  respectively;  the  right 
column  shows  the  10  mA  phases  after  4D-PWLS  pro¬ 
cessing.  From  top  to  bottom  are  phase  0%,  20%,  40%, 
60%,  80%,  respectively.  The  red  rectangles  represent 
the  selected  ROI  for  calculation  of  SNRs,  each  of  which 
contains  5X5X5  voxels. 
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Table  II.  A  comparison  of  SNRs  in  the  thorax  phantom  study  for  different  phases.  The  penalty  a  was  2.0 
X  10-6  in  the  4D-PWLS  smoothing. 


Phase  0% 

Phase  20% 

Phase  40% 

Phase  60% 

Phase  80% 

100  mA  scan 

0.479 

0.215 

0.122 

0.226 

0.349 

10  mA  scan 

0.078 

0.057 

0.015 

0.045 

0.059 

PWLS  Smoothed 

0.189 

0.170 

0.086 

0.167 

0.212 

10  mA  scan 

calculated  the  signal-to-noise  ratios  (SNRs),  defined  as 
MEAN/SD  (mean  divided  by  standard  deviation)  within  a 
selected  region  of  interest  shown  in  Fig.  8,  and  listed  them  in 
Table  II.  From  the  table,  it  is  seen  that  the  average  SNR  of 
10  mA  images  increased  by  more  than  three-fold  from  0.051 
to  0.165  after  the  proposed  post-processing.  The  average 
SNR  of  100  mA  phase  images  is  0.278. 

While  useful,  the  above  measure  of  SNR  is  not  sufficient 
to  characterize  the  goodness  of  a  smoothing  method.  To  as¬ 
sess  the  proposed  algorithm  accurately,  we  calculated  the 
relative  contrast  (RC)  in  addition  to  the  background  noise, 
similar  to  Ref.  19.  In  our  study,  the  RC  is  defined  as  the 
signal  difference  between  object  and  background.  Without 
loss  of  generality,  one  slice  location  of  the  three  datasets  at 
phase  0%  was  used  to  measure  the  RCs  and  background 
noises.  The  slices  are  shown  in  Fig.  9  for  the  three  datasets 
including  the  high-mA,  low-mA,  and  the  smoothed  images, 
from  top  to  bottom,  respectively.  The  object  signal  was  mea¬ 
sured  as  the  mean  CT  number  of  nine  selected  holes  indi¬ 
cated  by  the  yellow  rectangles  in  Fig.  9,  each  containing  4 
X  4  pixels;  the  background  signal  was  the  mean  CT  number 
of  a  selected  uniform  region  close  to  the  holes  containing 
20X20  pixels  (the  red  rectangle  in  Fig.  9).  The  RCs  and 
background  noise  were  measured  for  the  proposed  algorithm 
for  different  settings  of  the  smoothing  parameter  a ,  and  their 
relations  are  plotted  in  Fig.  10,  where  the  red  lines  represent 
values  for  the  100  mA  data,  and  the  diamond  marks  are  for 
the  smoothed  data  at  a  different  smoothing  level.  When  a 
=  0,  it  represents  the  original  low-mA  (10  mA)  data.  It  is 
found  that  the  background  noise  decreases  rapidly  with  a  at 
the  beginning  and  the  decreasing  rate  gets  slower  when  a 
increases  further,  while  the  contrast  is  reduced  almost  lin¬ 
early  with  a.  According  to  the  plots  shown  in  Fig.  10,  a  good 
choice  for  the  parameter  a  in  this  case  seems  to  be  between 

1  X  10“6  and  1  X  10  5  because  of  a  large  percentage  of  reduc¬ 
tion  of  the  noise  and  relatively  small  loss  of  the  contrast.  In 
Table  III,  we  list  the  measured  RC  and  noise  for  the  proposed 
4D  PWFS  method  with  an  “optimal”  a  empirically  chosen  as 

2  X  10-6.  Compared  with  the  original  low-mA  data,  the  back¬ 
ground  noise  was  reduced  approximately  2.7  times,  with 
only  1.3%  loss  of  the  contrast. 

B.  Patient  study 

The  patient  study  for  the  end-inspiration  phase  (phase 
0%)  and  end-expiration  phase  (phase  50%)  are  shown  in 
Figs.  11(a)  and  12(a),  respectively.  The  left  column  in  each 
figure  is  the  original  phase  image  obtained  from  GE  Advan¬ 


tage  Workstation,  and  the  right  column  is  the  image  after 
processing  with  the  proposed  4D-PWFS  enhancement.  From 
top  to  bottom  are  the  axial,  coronal,  and  sagittal  displays, 
respectively.  Successful  noise  suppression  is  observed  for 
both  phases  in  these  pictures.  With  of=1.0X  10-5  the  SNRs 
increased  from  2.204  to  4.558  for  the  end-expiration  phase, 
and  from  1.741  to  3.862  for  the  end-inspiration  phase  in  the 
selected  ROIs  (red  rectangles  in  the  figures,  each  consisting 
of  125  voxels).  In  Figs.  11(b)  and  12(b),  the  horizontal  pro¬ 
files  across  the  center  of  the  transaxial  slices  are  compared 
for  the  original  and  processed  images.  It  is  found  that  the 
4D-PWFS  processing  preserves  the  CT  values  and  edge  in¬ 
formation  very  well  while  reducing  the  noise.  This  patient 
study  was  performed  using  a  4D  protocol  using  x-ray  tube 
current  of  90  mA.  For  a  thorax  scan,  it  is  possible  that  the 
current  can  be  reduced.  However,  it  is  important  to  empha¬ 
size  that  the  proposed  technique  represents  an  independent 
way  to  reduce  the  patient  radiation  dose. 

In  order  to  demonstrate  the  role  of  the  temporal  informa¬ 
tion  in  suppressing  noise,  we  compared  the  same  algorithm 
with  and  without  the  aid  of  the  registered  phases.  By  assign¬ 
ing  the  weights  of  the  temporal  neighbors  to  zero,  a  conven¬ 
tional  3D-PWFS  algorithm  can  be  obtained,  which  smoothes 
an  image  based  only  on  the  3D  spatial  correlation.  Results  of 
the  end-expiration  phase  using  the  3D  PWFS  method  with 
same  a  are  shown  in  the  top  row  of  Fig.  13.  The  SNRs  for 
the  selected  ROI  increase  from  2.204  to  2.783,  compared  to 
4.558  of  the  4D-PWFS  method.  Furthermore,  for  a  qualita¬ 
tive  examination  of  the  spatial  resolution,  we  compared  3D- 
and  4D-PWFS  results  by  subtracting  from  each  of  them  the 
original  noisy  image.  The  resulting  difference  images  are 
shown  in  the  middle  and  bottom  rows  in  Fig.  13,  where 
better  edge  preservation  is  observed  for  the  4D-PWFS 
method  as  well. 

IV.  DISCUSSION  AND  CONCLUSION 

We  have  proposed  a  4D-PWLS  smoothing  method  to  re¬ 
duce  the  noise  in  4D-CT  images.  The  technique  allows  us  to 
obtain  high  quality  4D-CT  images  based  on  the  data  acquired 
with  a  low  x-ray  tube  current  at  individual  phases,  resulting 
in  a  significant  reduction  in  the  patient  radiation  dose. 
Through  deformable  registration,  the  method  incorporates 
the  information  of  different  phases  into  one  objective  func¬ 
tion  and  finds  the  optimal  estimation  of  the  true  image  in 
terms  of  the  least  square  metric  based  on  the  first-  and 
second-order  statistics  of  the  data.  In  this  method,  the  data  in 
the  temporal  domain  are  incorporated  into  the  penalty  term 
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Fig.  10.  The  relations  of  noise  and  relative  contrast  with  the  penalty  weight 
a.  The  noise  is  characterized  by  the  standard  deviation  of  the  CT  numbers  in 
the  uniform  region  shown  in  Fig.  9.  The  red  lines  in  the  figure  show  the 
noise  and  relative  contrast  of  the  100  mA  image.  The  blue  line  and  diamond 
marks  are  for  the  smoothed  images  with  different  a.  When  cr=0  it  represents 
the  10  mA  original  data. 


Fig.  9.  Axial  CT  images  of  the  thorax  phantom.  From  top  to  bottom  are  the 
100  mA,  10  mA,  and  4D-PWLS  smoothed  10  mA  images.  The  yellow 
squares  show  the  air  holes  for  measuring  the  “signal,”  the  red  square  shows 
the  region  for  measuring  the  “background.”  The  relative  contrast  is  defined 
as  the  CT  number  difference  between  the  signal  and  the  background. 


as  additional  neighbors,  and  their  “distances”  are  determined 
by  the  maximized  NCC  in  the  registration  step  to  accommo¬ 
date  possible  local  inaccuracies.  One  concern  of  this  method 
is  that  if  the  tumor  contour  is  compromised  after  the  defor¬ 
mation  because  of  the  existence  of  small  uncertainty  in  de¬ 
formable  registration.  The  experience  from  another  group20 
and  our  group21  have  indicated  that  generally  an  accuracy  of 
less  than  2-3  mm  is  achievable,  and  the  statistical  averaging 
strategy  as  described  above  should  further  reduce  the  uncer¬ 
tainty.  While  the  method  seems  to  be  very  robust,  better  reg¬ 
istration  will  definitely  improve  the  performance  of  the  pro¬ 
posed  technique. 


The  proposed  4D-PWLS  technique  has  been  validated  by 
phantom  experiments  and  applied  to  patient  data  as  well, 
where  it  is  observed  that  the  method  effectively  smoothes  out 
the  noise  and  leads  to  clear  improvements  of  the  image  qual¬ 
ity.  Basically,  the  method  integrates  ten  low-dose  phases  into 
one  phase,  hence  will  result  in  an  image  comparable  to  a 
ten-times  higher  dose  scan  in  an  ideal  situation  (perfect  reg¬ 
istration),  which  means  a  maximum  of  ten-fold  dose  could 
be  saved.  Note  that  the  4D  PWLS  is  essentially  a  low-pass 


Table  III.  A  comparison  of  relative  contrast  and  noise  in  the  thorax  phan¬ 
tom  study. 


100  mA  scan 

10  mA  scan 

Smoothed 

10  mA  scan 

Relative  contrast  (HU) 

793.60 

793.72 

783.95 

Background  noise  (HU) 

26.98 

133.48 

48.86 

Contrast  to  noise  ratio 

29.4 

5.9 

16.0 
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Fig.  11.  (a)  Patient  study  for  the  4D- 
PWLS  method  at  the  end-inspiration 
phase.  The  left  column  contains  the 
original  images  acquired  from  the  GE 
Advantage  Workstation,  and  the  right 
column  shows  the  image  after  4D- 
PWLS  processing.  The  red  rectangles 
represent  the  selected  ROI  for  the  cal¬ 
culation  of  SNRs,  each  of  which  con¬ 
tain  5X5X5  voxels,  (b).  A  compari¬ 
son  of  the  horizontal  profiles  across 
the  center  of  the  transaxial  images  in 
(a). 


filtering  method,  which  inevitably  reduces  the  image  spatial/ 
temporal  resolution  while  smoothing  the  noise.  The  tradeoff 
between  the  noise  and  resolution  can  be  controlled  by  select¬ 
ing  an  appropriate  value  for  the  penalty  weight  a,  so  that  a 
reasonably  enhanced  image  can  be  achieved,  as  demon¬ 


strated  in  this  work.  However,  the  way  to  determine  the  op¬ 
timal  value  for  the  parameter  a  is  yet  empirical.  The  overall 
evaluation  of  the  smoothing  method  eventually  depends  on 
the  performance  tests  of  tumor  detectability.  For  example,  to 
quantify  the  effectiveness  of  the  method  in  terms  of  the  clini- 
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Fig.  12.  (a)  Patient  study  for  the  4D- 
PWLS  method  at  the  end-expiration 
phase.  The  left  column  contains  the 
original  images  acquired  from  the  GE 
Advantage  Workstation,  and  the  right 
column  shows  the  image  after  4D- 
PWLS  processing.  The  red  rectangles 
represent  the  selected  ROI  for  calcula¬ 
tion  of  SNRs,  each  of  which  contain 
5X5X5  voxels,  (b).  A  comparison  of 
the  horizontal  profiles  across  the  cen¬ 
ter  of  the  transaxial  images  in  (a). 


cal  goal  in  practice,  receiver  operating  characteristic  (ROC) 

22 

studies  should  be  carried  out  with  large  amount  of  data, 
which  is  beyond  the  scope  of  the  current  work.  Our  main 
purpose  in  this  paper  is  to  present  a  statistic  framework  to 
incorporate  temporal  information  into  image  restoration.  The 


ICM  algorithm  in  our  studies  usually  converged  to  a  satisfied 
solution  within  ten  iterations.  However,  because  of  a  qua¬ 
dratic  prior  being  used,  other  algorithms  such  as  a  conjugate 
gradient  may  be  developed  for  even  better  convergence.  Fur¬ 
thermore,  it  may  be  of  research  interest  to  develop  a  penalty 
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Fig.  13.  A  comparison  of  PWLS  with 
and  without  registered  phase  informa¬ 
tion.  The  top  row  is  the  result  of  3D 
PWLS  smoothing  without  using  tem¬ 
poral  information  for  the  end- 
inspiration  phase  patient  data.  The 
middle  row  shows  the  difference  im¬ 
ages  between  the  3D-PWLS  smooth¬ 
ing  results  (top  row  in  this  figure)  and 
original  images  [the  left  column  in 
Fig.  11(a)].  The  bottom  row  contains 
the  difference  images  between  4D- 
PWLS  smoothing  results  [the  right 
column  in  Fig.  11(a)]  and  original  im¬ 
ages.  The  fewer  edges  observed  in  the 
difference  image  indicate  an  improved 
preservation  of  the  spatial  resolution 
in  the  4D  method.  The  red  rectangles 
represent  the  selected  ROI  for  calcula¬ 
tion  of  SNRs,  each  of  which  contain 
5X5X5  voxels. 


term  other  than  the  quadratic  form  for  better  preservation  of 
the  spatial/temporal  resolution  for  4D  CT. 

In  our  earlier  study,15  it  was  found  that  further  improve¬ 
ments  to  image  quality  may  be  possible  if  the  noise  reduction 
is  performed  in  sinogram  space  (before  reconstruction)  rather 
than  in  image  space  (after  reconstruction).  In  future  research, 
we  will  study  the  possibility  of  applying  deformable  regis¬ 
tration  and  4D-PWLS  to  sinogram  space.  A  quantitative 
evaluation  of  the  benefits  of  the  proposed  method  to  appli¬ 
cations  in  PET/CT  imaging23-25  and  4D  treatment  planning 
for  radiotherapy26,27  are  also  under  investigation. 

ACKNOWLEDGMENTS 

This  research  is  supported  in  part  by  the  Department  of 
Defense  (DAMD17-03-1-0657)  and  the  National  Cancer  In¬ 
stitute  (5  ROI  CA98523-01). 

a) Author  to  whom  correspondence  should  be  addressed.  Department  of 
Radiation  Oncology,  Stanford  University  School  of  Medicine,  Clinical 
Cancer  Center,  875  Blake  Wilbur  Drive,  Rm  CC-G204,  Stanford,  Cali¬ 
fornia  94305-5847.  Telephone:  (650)  498-7896;  Fax:  (650)  498-4015; 
Electronic  mail:  lei@reyes.stanford.edu 
‘s.  S.  Vedam,  P.  J.  Keall,  V.  R.  Kini,  H.  Mosafavi,  H.  P.  Shukla,  and  R. 
Mohan,  “Acquiring  a  four-dimensional  computed  tomography  dataset  us¬ 
ing  an  external  respiratory  signal,”  Phys.  Med.  Biol.  48,  45-62  (2003). 
2P.  J.  Keall,  G.  Starkschall,  H.  Shukla,  K.  M.  Forster,  V.  Ortiz,  C.  W. 
Stevens,  S.  S.  Vedam,  R.  George,  T.  Guerrero,  and  R.  Mohan,  “Acquiring 
4D  thoracic  CT  scans  using  a  multislice  helical  method,”  Phys.  Med. 
Biol.  49,  2053-2067  (2004). 

3D.  A.  Low,  M.  Nystrom,  E.  Kalinin  et  al.,  “A  method  for  the  reconstruc¬ 
tion  of  four-dimensional  synchronized  CT  scans  acquired  during  free 
breathing,”  Med.  Phys.  30,  1254-1263  (2003). 

4T.  Pan,  T.-Y.  Lee,  E.  Rietzel,  and  G.  T.  Y.  Chen,  “4D-CT  imaging  of  a 
volume  influenced  by  respiratory  motion  on  multi-slice  CT,”  Med.  Phys. 
31,  333-340  (2004). 


5Y.  E.  Erdi,  S.  A.  Nehmeh,  T.  Pan,  S.  Kohlmyer,  O.  D.  Squire,  H.  Schoder, 
H.  W.  Yeung,  J.  L.  Humm,  K.  E.  Rosenzweig,  and  S.  M.  Larson,  “The  CT 
motion  quantitation  of  lung  lesions  and  its  impact  on  PET  measured 
SUV’s, ”J.  Nucl.  Med.  45,  1287-1292  (2004). 

6J.  Lian,  L.  Xing,  S.  Hunjan,  C.  Dumoulin,  J.  Levin,  A.  Lo,  R.  Watkins,  K. 
Rohling,  R.  Giaquinto,  D.  Kim,  D.  Spielman,  and  B.  Daniel,  “Mapping  of 
the  prostate  in  endorectal  coil-based  MRI/MRSI  and  CT:  A  deformable 
registration  and  validation  study,”  Med.  Phys.  31,  3087-3094  (2004). 

7E.  Schreibmann  and  L.  Xing,  “Narrow  band  deformable  registration  of 
prostate  MRI/MRSI  and  CT  studies,”  Int.  J.  Radiat.  Oncol.,  Biol.,  Phys., 
62,  595-605  (2005). 

8D.  Mattes,  D.  R.  Haynor,  H.  Vesselle,  T.  K.  Lewellen,  and  W.  Eubank, 
“PET-CT  image  registration  in  the  chest  using  free-form  deformations,” 
IEEE  Trans.  Med.  Imaging  22,  120-128  (2003). 

9T.  Rohlfing,  C.  R.  Maurer,  Jr.,  D.  A.  Bluemke,  and  M.  A.  Jacobs, 
“Volume-preserving  nonrigid  registration  of  MR  breast  images  using  free¬ 
form  deformation  with  an  incompressibility  constraint,”  IEEE  Trans. 
Med.  Imaging  22,  730-741  (2003). 

10D.  Rueckert,  L.  I.  Sonoda,  C.  Hayes,  D.  L.  Hill,  M.  O.  Leach,  and  D.  J. 
Hawkes,  “Nonrigid  registration  using  free-form  deformations:  Applica¬ 
tion  to  breast  MR  images,”  IEEE  Trans.  Med.  Imaging  18,  712-721 
(1999). 

nD.  Loeckx,  F.  Maes,  D.  Vandermeulen,  and  P.  Suetens,  “Non-rigid  image 
registration  using  a  statistical  spline  deformation  model,”  Info.  Process. 
Med.  Imaging  18,  463-474  (2003). 

12D.  C.  Liu  and  J.  Nocedal,  “On  the  limited  memory  BFGS  method  for 
large  scale  optimization,”  Math.  Program.  45,  503-528  (1989). 

13B.  F.  Hutton  and  M.  Braun,  “Software  for  image  registration:  Algorithm, 
accuracy,  efficacy,”  Semin.  Nucl.  Med.  33,  180-192  (2003). 

14J.  A.  Fessler,  “Penalized  weighted  least-squares  image  reconstruction  for 
positron  emission  tomography,”  IEEE  Trans.  Med.  Imaging  13,  290-300 
(1994). 

15T.  Li,  X.  Li,  J.  Wang,  J.  Wen,  H.  Lu,  J.  Hsieh,  and  Z.  Liang,  “Nonlinear 
sinogram  smoothing  for  low-dose  x-ray  CT,”  IEEE  Trans.  Nucl.  Sci.  51, 
2505—2513  (2004). 

16S.  Geman  and  D.  Geman,  “Stochastic  relaxation,  Gibbs  distributions  and 
the  Bayesian  restoration  of  images,”  IEEE  Trans.  Pattern  Anal.  Mach. 
Intell.  6,  721-741  (1984). 

17X.  Li,  L.  Li,  H.  Lu,  and  Z.  Liang,  “Partial  volume  segmentation  of  brain 
magnetic  resonance  images  based  on  maximum  a  posteriori  probability,” 


Medical  Physics,  Vol.  32,  No.  12,  December  2005 


3660 


Li  et  a!.:  Radiation  dose  reduction  in  4D  CT 


3660 


Med.  Phys.  32,  2337-2345  (2005). 

18J.  Besag,  “On  the  statistical  analysis  of  dirty  pictures,”  J.  R.  Stat.  Soc.  Ser. 
B.  Methodol.  48,  259-302  (1986). 

19P.  K.  Saha  and  J.  K.  Udupa,  “Scale-based  diffusive  image  filtering  pre¬ 
serving  boundary  sharpness  and  fine  structures,”  IEEE  Trans.  Med.  Im¬ 
aging  20,  1140-1155  (2001). 

20M.  M.  Coselmon,  J.  M.  Balter,  D.  L.  McShan,  and  M.  L.  Kessler,  “Mu¬ 
tual  information  based  CT  registration  of  the  lung  at  exhale  and  inhale 
breathing  states  using  thin-plate  splines,”  Med.  Phys.  31,  2942-2948 
(2004). 

21E.  Schreibmann,  G.  T.  Y.  Chen,  and  L.  Xing,  “Image  interpolation  in  4D 
CT  using  a  BSpline  deformable  registration  model,”  Int.  J.  Radiat.  Oncol., 
Biol.,  Phys.. 

22T.  Li,  J.  Wen,  G.  Han,  H.  Lu,  and  Z.  Liang,  “Evaluation  of  an  efficient 
compensation  method  for  quantitative  fan-beam  brain  SPECT  reconstruc¬ 
tion,”  IEEE  Trans.  Med.  Imaging  24,  170-179  (2005). 

23S.  A.  Nehmeh,  Y.  E.  Erdi,  T.  Pan,  A.  Pevsner,  K.  E.  Rosenzweig,  E. 
Yorke,  G.  S.  Mageras,  H.  Schoder,  Phil  Vernon,  O.  Squire,  H.  Mostafavi, 


S.  M.  Larson,  and  J.  L.  Humm,  “Four-dimensional  (4D)  PET/CT  imaging 
of  the  thorax,”  Med.  Phys.  31,  3179-3186  (2004). 

24S.  A.  Nehmeh,  Y.  E.  Erdi,  K.  E.  Rosenzweig,  H.  Schoder,  S.  M.  Larson, 
O.  D.  Squire,  and  J.  L.  Humm,  “Reduction  of  respiratory  motion  artifacts 
in  PET  imaging  of  lung  cancer  by  respiratory  correlated  dynamic  PET: 
Methodology  and  comparison  with  respiratory  gated  PET,”  J.  Nucl.  Med. 
44,  1644-1648  (2003). 

25S.  A.  Nehmeh,  Y.  E.  Erdi,  T.  Pan,  E.  Yorke,  G.  S.  Mageras,  K.  E.  Rosen¬ 
zweig,  H.  Schoder,  H.  Mostafavi,  O.  D.  Squire,  A.  Pevsner,  S.  M.  Larson, 
and  J.  L.  Humm,  “Quantitation  of  respiratory  motion  during  4D-PET/CT 
acquisition,”  Med.  Phys.  31,  1333-1338  (2004). 

26P.  J.  Keall,  S.  Vedam,  V.  Kini,  S.  Joshi,  G.  Tracton,  and  R.  Mohan,  “4D 
IMRT:  Planning  methodology,”  Int.  J.  Radiat.  Oncol.,  Biol.,  Phys.  54, 
S318-S319  (2002). 

27H.  Shirato,  S.  Shimizu,  K.  Kitamura  et  al.,  “Four-dimensional  treatment 
planning  and  fluoroscopic  real-time  tumor  tracking  radiotherapy  for  mov¬ 
ing  tumor,”  Int.  J.  Radiat.  Oncol.,  Biol.,  Phys.  48,  435-442  (2000). 


Medical  Physics,  Vol.  32,  No.  12,  December  2005 


Institute  of  Physics  Publishing 


Physics  in  Medicine  and  Biology 


Phys.  Med.  Biol.  50  (2005)  1469-1482 


doi :  1 0. 1 088/003 1-915 5/5 0/7/0 1 0 


Quantitation  of  the  a  priori  dosimetric  capabilities  of 
spatial  points  in  inverse  planning  and  its  significant 
implication  in  defining  IMRT  solution  space* 


Z  Shou1,  Y  Yang1,  C  Cotrutz1,  D  Levy2  and  Lei  Xing1 

1  Department  of  Radiation  Oncology,  Stanford  University,  Stanford,  CA  94305-5847,  USA 

2  Department  of  Mathematics,  Stanford  University,  Stanford,  CA  94305-2125,  USA 

E-mail:  lei@reyes.stanford.edu 

Received  4  January  2005,  in  final  form  2  February  2005 

Published  16  March  2005 

Online  at  stacks. iop.org/PMB/50/ 1469 

Abstract 

In  inverse  planning,  the  likelihood  for  the  points  in  a  target  or  sensitive  structure 
to  meet  their  dosimetric  goals  is  generally  heterogeneous  and  represents  the 
a  priori  knowledge  of  the  system  once  the  patient  and  beam  configuration  are 
chosen.  Because  of  this  intrinsic  heterogeneity,  in  some  extreme  cases,  a  region 
in  a  target  may  never  meet  the  prescribed  dose  without  seriously  deteriorating 
the  doses  in  other  areas.  Conversely,  the  prescription  in  a  region  may  be 
easily  met  without  violating  the  tolerance  of  any  sensitive  structure.  In  this 
work,  we  introduce  the  concept  of  dosimetric  capability  to  quantify  the  a  priori 
information  and  develop  a  strategy  to  integrate  the  data  into  the  inverse  planning 
process.  An  iterative  algorithm  is  implemented  to  numerically  compute  the 
capability  distribution  on  a  case  specific  basis.  A  method  of  incorporating  the 
capability  data  into  inverse  planning  is  developed  by  heuristically  modulating 
the  importance  of  the  individual  voxels  according  to  the  a  priori  capability 
distribution.  The  formalism  is  applied  to  a  few  specific  examples  to  illustrate  the 
technical  details  of  the  new  inverse  planning  technique.  Our  study  indicates  that 
the  dosimetric  capability  is  a  useful  concept  to  better  understand  the  complex 
inverse  planning  problem  and  an  effective  use  of  the  information  allows  us  to 
construct  a  clinically  more  meaningful  objective  function  to  improve  IMRT 
dose  optimization  techniques. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


*  Part  of  this  work  was  presented  in  the  14th  International  Conference  on  the  Use  of  Computers  in  Radiation  Therapy, 
Seoul,  Korea,  2004. 


0031-9155/05/071469+14$30.00  ©  2005  IOP  Publishing  Ltd  Printed  in  the  UK 


1469 


1470 


Z  Shou  et  al 


1.  Introduction 

One  of  the  implicit  assumptions  in  current  inverse  planning  is  that  all  points  within 
a  target  or  sensitive  structure  are  equivalent  (Brahme  et  al  1982,  Bortfeld  et  al  1990, 
Cho  et  al  1998,  Gopal  and  Starkschall  2001,  Holmes  and  Mackie  1994,  Langer  and  Leong 
1987,  Spirou  and  Chui  1998,  Webb  1991,  Xing  and  Chen  1996,  Zagars  et  al  2002).  In  reality, 
not  all  voxels  have  the  same  chance  of  complying  with  the  prescription  because  the  dose- 
limiting  factors  imposed  by  the  involved  sensitive  structures  are  not  uniformly  distributed  in 
space  and  the  dose  delivery  is  depth  dependent.  For  a  given  patient  and  beam  configuration, 
it  is  practically  useful  to  know  the  likelihood  for  each  individual  point  within  a  target  or 
sensitive  structure  to  meet  its  dosimetric  goal.  By  understanding  this  intrinsic  property  of 
the  system,  one  can  better  model  the  therapeutic  plan  optimization  problem  and  improve  the 
inverse  planning  techniques. 

In  this  work,  we  introduce  the  concept  of  dosimetric  capability  for  an  arbitrary  voxel  in  a 
target  or  sensitive  structure  to  quantify  the  likelihood  for  the  voxel  to  meet  the  specified  dose. 
The  capability  calculation  finds  the  potentially  problematic  regions  and,  more  importantly, 
the  degree  of  problems  for  these  regions  to  meet  their  goals,  and  permits  us  to  purposely 
modify  the  penalty  strategy  during  the  construction  of  objective  function  to  minimize  the 
problem.  Mathematically,  this  process  is  realized  by  heuristically  modulating  the  importance 
of  the  individual  voxels  according  to  the  a  priori  capability  distribution.  The  inverse  planning 
formalism  with  dosimetric  capability-modulated  importance  factors  is  applied  to  a  few  specific 
examples  to  illustrate  the  technical  details.  The  approach  sheds  useful  insight  into  the 
inverse  planning  problem  and  allows  us  to  search  for  IMRT  solutions  that  would  otherwise  be 
inaccessible. 


2.  Methods  and  materials 


2.1.  Definition  of  dosimetric  capability 


We  quantify  the  likelihood  for  a  voxel  to  meet  its  dosimetric  goal  by  introducing  the  concept  of 
dosimetric  capability.  Let  us  first  consider  the  case  of  one  incident  beam.  For  a  target  voxel,  the 
capability  is  assessed  by  the  degree  to  which  the  voxel  meets  the  prescription  without  violating 
the  tolerance  of  the  sensitive  structure.  The  maximum  achievable  dose,  Dach(ncr),  at  the  voxel 
na  in  the  target  is  determined  by  scaling  the  intensity  of  the  contributing  beamlet  to  the  highest 
value  set  by  the  tolerance  of  the  sensitive  structure.  Mathematically,  the  ‘capability’,  r /,  is 
defined  as 


r]{na) 


D*c\na) 

Dr 


(1) 


The  evaluation  of  equation  (1)  is  straightforward  for  the  case  of  a  single  incident  beam 
or  when  there  is  no  overlap  of  beamlets  at  the  dose-limiting  voxel  in  the  sensitive  structure. 
For  multi-field  IMRT,  the  Ddch(nc r)  is  determined  not  only  by  those  beamlets  that  directly 
intercept  the  voxel  n ,  but  also  possibly  by  other  beamlets  irradiating  different  parts  of  the 
target.  The  coupled  system  can  be  described  by  a  set  of  linear  equations  (which  can  be  easily 
written  from  the  above  definition  of  the  maximum  achievable  dose  and  the  dose  as  a  function 
of  beamlet  weights)  and  the  Dach(na)  of  a  target  voxel  will  be  obtained  by  optimizing  the 
linear  system,  under  the  condition  that  the  dose  at  any  sensitive  structure  voxel  is  equal  to 
its  tolerance.  The  beam  profiles  so  obtained  deliver  the  maximum  dose  to  the  target  without 
violating  the  tolerances  of  the  sensitive  structures.  Mathematically,  the  system  equations  are 
underdetermined  and  a  Cimmino  algorithm  (Stark  and  Yang  1988,  Xiao  et  al  2000),  which 
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was  first  applied  to  radiation  therapy  by  Starkschal  and  Eifel  (1992),  can  be  used  to  find  the 
solution.  After  assigning  all  beamlets  with  high  enough  initial  intensities  (say,  two  or  three 
times  the  intensity  values  that  deliver  the  prescribed  dose  to  the  target)  and  dose  calculation, 
the  calculation  consists  of  the  following  steps. 


(a)  Choose  a  beamlet  in  a  beam  and  locate  the  voxels  in  sensitive  structures  that  are  traversed 
by  the  beamlet. 

(b)  For  each  located  voxel  in  the  sensitive  structures,  check  if  the  tolerance  is  exceeded.  If 
yes,  decrease  the  value  of  the  beamlet  based  on 


with 

Cn„  (U>k) 


Wjm  ~  W  jm  +  ^  ^  '  Cn„  (w  )djm(n< r)> 
na 


(2) 


(£>r  -  Pg(nff)) 

EjmdjjHo) 

0 


if  Dkc(na)  >  D^x  and  na  E  sensitive  structures 
otherwise. 


(c)  Update  the  dose. 

(d)  Repeat  steps  (b)  and  (c)  until  all  voxels  intercepted  by  the  beamlet  are  checked. 

(e)  Repeat  (a)-(d)  for  the  next  beamlet. 


The  above  calculation  is  repeated  until  the  doses  in  sensitive  structures  are  equal  to  their 
tolerances.  The  relaxation  parameter  X  is  generally  set  to  a  small  value  (0.02  <  X  <  0.1)  to 
ensure  a  smooth  updating  of  the  beamlet  weights.  The  beamlet-by-beamlet  updating  scheme 
is  similar  to  the  simultaneous  iterative  inverse  planning  technique  (SIITP)  (Xing  and  Chen 
1996).  However,  the  goal  here  is  to  search  for  the  beam  profiles  that  deliver  the  highest 
achievable  doses  in  the  target,  Z)ach  (na),  without  violating  the  dose  tolerance  of  the  involved 
sensitive  structures  (in  other  words,  any  increase  in  the  beamlet  weights  would  lead  to  a  dose 
exceeding  the  tolerance  in  one  or  multiple  points  inside  a  sensitive  structure).  Mathematically, 
the  above  calculation  is  an  iterative  projection  algorithm  converging  to  a  feasible  solution  if 
the  constraints  are  satisfied,  otherwise  to  a  compromise  solution  that  minimizes  the  weighted 
sum  of  squares  of  deviations  from  the  tolerance  doses  if  dose  constraints  are  violated. 

The  dosimetric  capability  of  a  voxel  in  a  sensitive  structure  is  characterized  by 
the  minimum  achievable  dose  at  the  voxel  needed  in  order  to  administer  the  prescribed  dose 
to  the  target  voxels.  For  consistency,  we  denote  the  minimum  achievable  dose  by  D*ch(na) 
(where  na  represents  a  voxel  in  the  sensitive  structure)  and  quantify  the  dosimetric  capability 
of  the  voxel  according  to 


rj(na) 


Dach(na) ' 


(3) 


For  a  voxel  in  a  sensitive  structure,  the  lower  the  minimum  achievable  dose,  the  more 
‘capable’  the  voxel  is.  The  evaluation  of  equation  (3)  is  once  again  straightforward  for  the 
case  with  a  single  incident  beam.  To  obtain  the  minimum  achievable  dose,  we  first  set  the 
corresponding  beamlet  to  such  an  intensity  that  the  prescribed  dose  is  delivered  to  the  target 
voxels,  and  then  evaluate  the  dose,  Dach(na),  at  na  in  the  sensitive  structure.  For  the  case  of 
multiple  incident  beams,  a  set  of  linear  equations  needs  to  be  solved  under  the  condition  that  the 
target  voxels  receive  their  prescription  doses.  Similar  to  the  target  case,  the  system  equations 
become  undetermined  when  there  are  multiple  incident  beams  and  a  Cimmino  algorithm  is 
used  to  find  the  beam  profiles  that  yield  the  lowest  achievable  dose  in  the  sensitive  structures. 
We  first  assign  all  beamlets  with  low  enough  initial  intensities  so  that  doses  delivered  to  the 
target  voxels  are  less  than  the  prescription,  and  then  do  the  following. 


1472 


Z  Shou  et  al 


(a)  Choose  a  beamlet  in  a  beam  and  locate  the  voxels  in  the  target  that  are  traversed  by  the 
beamlet. 

(b)  For  each  located  voxel  in  the  target,  check  if  the  dose  prescription  is  exceeded.  If  yes, 
lower  the  beamlet  weight  based  on  equation  (2)  with 


C„0(vuk ) 


(£>r  -  P*(ng)) 
Ejmdjmina) 

0 


if  Dkc(na)  <  D„k  and  na  e  target 
otherwise. 


(c)  Update  the  doses. 

(d)  Repeat  steps  (b)  and  (d)  until  all  voxels  intercepted  by  the  beamlet  are  checked. 

(e)  Repeat  (a)-(d)  for  the  next  beamlet. 

The  above  calculation  proceeds  iteratively  until  the  target  dose  prescription  is  completely 
met  for  all  voxels.  We  emphasize  that  the  beam  profiles  obtained  above  are  not  intended  to 
approximate  the  optimal  fluence  profiles  for  IMRT  treatment.  Instead,  they  are  obtained  purely 
for  the  purpose  of  evaluating  the  dosimetric  capabilities  of  voxels  in  a  target  or  a  sensitive 
structure.  Once  the  capability  maps  are  obtained,  they  can  be  incorporated  into  an  inverse 
planning  procedure. 


2.2.  Incorporating  dosimetric  capability  into  inverse  planning 


The  capability  distribution  contains  information  about  the  intrinsic  heterogeneity  of  the 
dosimetric  capability  and  reveals  which  points  are  likely  to  violate  the  prescription.  The 
information  provides  a  guiding  map  for  the  inverse  planning  algorithm  to  differentially  deal 
with  the  regions  having  different  chances  of  meeting  their  dosimetric  goals.  Our  strategy  is 
to  assign  a  higher  penalty  (or  high  local  importance)  to  those  voxels  with  lower  capabilities 
(those  voxels  are  likely  to  have  a  dose  lower  than  the  prescription  if  they  are  in  the  target 
volume,  or  higher  than  the  tolerance  if  they  are  located  in  a  sensitive  structure).  While  it  is  not 
difficult  to  intuitively  conceive  the  general  behaviour  of  the  relation  between  local  importance 
and  the  capability,  the  specific  form  of  the  relation  is  a  matter  of  experimenting.  The  relation 
between  the  local  importance  and  the  capability  used  in  our  study  is 

r(na)  =  (1  +g(na)2)ra,  (4) 


where  g(na)  is  empirically  defined  as 
3 


g(na)  = 


(cax  - 1) 

3 

(cm  - ') 


(j]ina)  -  1) 


(n(na)  -  1) 


na  e  target  and  rj(na)  ^  1 

na  g  target  or  sensitive  structure  and  rj{na)  <  1. 

(5) 


The  capability  map  and  the  corresponding  ^ax  or  are  obtained  for  each  structure. 
For  a  sensitive  structure,  if  r]{na)  >  1,  it  means  that  the  voxel  has  no  difficulty  meeting  its 
dosimetric  goal  and  g(na)  is  set  to  zero.  The  differential  penalty  scheme  allows  the  system 
to  suppress  potential  hot  spots  in  the  sensitive  structures  and  boost  the  potential  cold  spots  in 
the  target  volume  so  that  a  more  uniform  dose  distribution  can  be  achieved  in  the  target  while 
sparing  more  sensitive  structures. 
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Figure  1.  Dosimetric  capability  and  importance  maps  of  the  target  and  sensitive  structure  for  a 
hypothetical  case  for  two  beam  configurations:  (a)  the  single  beam  and  (b)  five  equally  spaced 
beams.  The  data  for  each  structure  are  normalized  to  unity.  For  visual  purposes,  the  capability 
and  importance  maps  of  the  target  and  sensitive  structure  are  enlarged  and  shown  in  the  left  and 
right  panels  for  each  beam  configuration.  The  lower  middle  panel  of  each  set  of  figures  shows  the 
complete  geometry  of  the  hypothetical  structures. 


2.3.  Implementation 

We  implemented  a  software  module  to  optimize  the  system  in  the  platform  of  the  PLUNC 
treatment  planning  system  (an  open  source  treatment  planning  system  from  the  University  of 
North  Carolina,  Chapel  Hill,  NC).  The  dose  calculation  engine  and  varieties  of  plan  evaluation 
tools  of  the  PLUNC  system  are  used  to  evaluate  and  compare  the  optimization  results.  The 
SIITP  (Xing  and  Chen  1996)  was  employed  to  obtain  the  optimal  beam  intensity  profiles.  The 
final  IMRT  plan  is  obtained  in  a  similar  manner  to  conventional  inverse  planning  except  that 
the  uniform  importance  for  the  involved  structures  is  replaced  by  the  non-uniform  importance 
distributions  given  by  equations  (4)  and  (5).  The  calculation  was  performed  on  a  PC  with  P4 
1.7  GHz  and  1024  MB  RAM. 

2.4.  Test  of  the  new  dose  optimization  formalism 

To  better  understand  the  physics  behind  the  capability  calculation,  we  first  constructed  a 
hypothetical  surrogate  case  (figure  1)  and  studied  the  behaviour  of  the  system  using  a  single 
beam  and  five  incident  beams.  The  gantry  angle  used  for  target  irradiation  in  the  single 
beam  case  was  320°,  and  in  the  five-beam  case  the  angles  used  were  32°,  104°,  176°,  248° 
and  320°,  where  IEC  convention  for  gantry  angle  is  used.  The  incident  photon  energy  was 
15  MV.  In  both  cases,  the  target  was  prescribed  to  100  (arbitrary  units)  and  the  sensitive 
structure  tolerance  was  set  to  be  25.  An  IMRT  treatment  was  also  planned  with  the  same 
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Table  1.  Summary  of  the  parameters  used  for  planning  the  IMRT  prostate  treatment.  The 
tolerances  of  the  sensitive  structures  are  used  in  the  evaluation  of  the  capability  maps. 


Prostate 

Bladder 

Rectum 

Femoral  heads 

Skin 

Importance  factor 

0.2 

0.05 

0.1 

0.05 

0.6 

Prescription/tolerance 

100 

65 

60 

70 

85 

Table  2.  Summary  of  the  parameters  used  for  planning  the  IMRT  treatment  of  the  paraspinal 
tumour.  The  tolerances  of  the  involved  sensitive  structures  are  used  in  the  evaluation  of  the 
capability  maps. 


GTV 

Spinal  cord 

Liver 

Kidney 

Tissue 

Importance  factor 

0.86 

0.03 

0.005 

0.05 

0.055 

Prescription/tolerance 

100 

30 

40 

30 

75 

five-beam  configuration  but  structurally  uniform  importance  factors,  and  the  result  was 
compared  with  the  newly  obtained  plan.  To  ensure  a  fair  plan  comparison,  in  this  and 
following  examples  the  importance  factors  were  chosen  in  such  a  way  that  the  target  DVHs 
were  the  same  for  the  cases  with  uniform  and  non-uniform  importance  factors.  The  net 
improvements  can  then  be  assessed  by  the  doses  given  to  the  sensitive  structures. 

The  new  algorithm  was  also  applied  to  two  clinical  cases:  a  prostate  case  and  a  paraspinal 
tumour  treatment.  To  illustrate  the  advantage  of  the  technique,  the  results  were  compared  with 
those  obtained  using  the  conventional  inverse  planning  with  uniform  importance  factors.  For 
the  prostate  IMRT  case,  six  equally  spaced  beams  starting  from  0°  were  used.  Some  relevant 
parameters  used  for  planning  the  patient  are  summarized  in  table  1.  As  in  conventional 
inverse  planning,  the  structure- specific  importance  factors,  {ra},  were  determined  by  trial- 
and-error.  To  assess  the  new  technique,  we  planned  the  case  using  three  different  strategies: 
(i)  uniform  importance  for  the  prostate  target  and  the  sensitive  structures;  (ii)  dosimetric 
capability-based  non-uniform  importance  for  the  prostate  target  and  uniform  importance  for 
all  sensitive  structures  and  (iii)  dosimetric  capability-based  non-uniform  importance  for  both 
prostate  target  and  the  sensitive  structures.  The  three  plans  were  compared  using  DVHs  and 
isodose  distribution  plots. 

In  the  paraspinal  case,  the  sensitive  structures  involved  include  the  spinal  cord,  liver 
and  kidney.  In  this  study,  five  15  MV  non-equally  spaced  coplanar  beams  (40°,  110°,  180°, 
255°  and  325°)  were  used  for  the  treatment.  The  structure- specific  importance  factors  and 
prescription/ tolerances  are  summarized  in  table  2.  IMRT  plans  with  and  without  importance 
factor  modulation  were  obtained  and  quantitative  comparison  was  performed. 

3.  Results  and  discussions 

3.1.  A  hypothetical  test  case 

The  capability  and  importance  distributions  for  both  single  beam  and  five-beam  configurations 
are  plotted  in  figure  1 .  The  isocapability  and  iso -importance  curves  are  normalized  to  unity. 
For  the  single  beam  calculation  (figure  1(a)),  it  is  seen  from  the  capability  map  (left  panel) 
that  the  target  region  in  the  middle  of  the  incident  beam  has  lower  dosimetric  capability  due 
to  the  restriction  of  the  tolerance  of  the  sensitive  structure.  For  the  points  on  both  sides  of 
the  central  target  region,  the  dosimetric  capabilities  are  much  higher  because  of  the  absence 
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Figure  2.  Isodose  distributions  for  plans  obtained  without  (left)  and  with  (right)  local  importance 
factor  modulation.  The  relative  isodose  curves  labelled  in  the  plots  are,  from  the  centre,  105% 
(red),  100%  (pink),  80%  (yellow)  and  40%  (blue),  respectively. 


of  sensitive  structure  ‘blocking’.  The  capability  distribution  within  the  sensitive  structure  can 
also  be  intuitively  interpreted.  For  the  voxels  distant  from  the  target  volume,  the  capability  is 
relatively  high,  indicating  that  these  voxels  are  less  dose-limiting  points  in  comparison  with 
the  voxels  close  to  the  target.  For  the  single  beam  case,  the  importance  map  (right  panel  of 
figure  1(a))  is  almost  an  inversion  of  the  capability  map. 

Figure  1(b)  shows  the  capability  and  importance  distributions  when  five  incident  beams 
are  used  to  irradiate  the  target.  The  calculation  is  fairly  efficient;  it  took  less  than  6  min  for 
the  system  to  obtain  the  capability  and  importance  maps.  In  this  case,  the  low  capability 
region  in  the  middle  of  the  target  shown  in  figure  1(a)  disappears  and  only  a  few  isolated  low 
capability  spots  show  up  near  the  edge  of  the  target.  On  the  other  hand,  the  overall  behaviour 
of  the  capability  map  in  the  sensitive  structure  is  not  changed  dramatically.  The  dosimetric 
capabilities  for  those  points  close  to  the  target  remain  relatively  low,  which  is  consistent  with 
our  intuition  since,  relatively  speaking,  these  points  are  more  dose-limiting  compared  to  the 
points  far  away  from  the  target.  The  results  also  suggest  that  the  boundary  region  between  the 
two  structures  is  likely  to  be  underdosed  (for  target)  or  overdosed  (for  sensitive  structure).  In 
order  to  improve  this,  a  non-uniform  penalty  scheme  derived  a  priori  or  a  posteriori  becomes 
necessary.  The  technique  described  in  this  paper  represents  an  example  of  the  a  priori 
method,  whereas  an  adaptive  adjustment  of  the  local  importance  factor  distribution  would  be 
an  example  of  the  latter.  In  general,  the  importance  distribution  is  not  a  simple  inversion  of 
the  capability  map  and  depends  also  on  the  absolute  values  of  the  capabilities,  as  suggested 
by  equations  (4)  and  (5).  This  is  especially  true  for  the  target,  in  which  we  need  not  only  to 
‘boost’  the  potential  cold  region(s)  but  also  to  ‘suppress’  the  potential  hot  spot(s).  The  high 
importance  in  the  left-middle  region  of  the  target  (see  the  right  panel  of  figure  1(b))  is  a  direct 
consequence  of  the  stated  requirement.  Indeed,  a  careful  examination  of  the  target  capability 
data  indicated  that  the  region  has  high  dosimetric  capability  and  is  thus  likely  to  be  overdosed. 
The  assignment  of  higher  local  importance  according  to  equation  (4)  permits  us  to  suppress 
the  potential  overdosing  a  priori. 

In  figure  2,  we  show  the  IMRT  plans  obtained  without  and  with  modulating  the  spatial 
importance  distribution.  The  left  panel  is  the  conventional  IMRT  plan  with  uniform  importance 
factors  assigned  to  both  target  and  sensitive  structures.  The  right  panel  shows  the  plan  obtained 
with  the  spatially  non-uniform  importance  distribution  plotted  in  the  right  panel  of  figure  1(b). 
It  is  clearly  seen  that  the  isodose  curves  in  the  sensitive  structure  are  ‘pushed’  towards  the  target 
and  the  dose  gradient  at  the  tumour  boundary  is  greatly  increased.  The  significant  improvement 
can  also  be  seen  in  the  DVH  plot  (figure  3).  It  is  remarkable  that,  by  simply  modulating  the 
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Figure  3.  Target  and  sensitive  structure  DVHs  corresponding  to  plans  with  two  different  penalty 
schemes:  (a)  uniform  importance  for  every  structure  (dotted-dash  curves)  and  (b)  non-uniform 
importance  for  both  target  and  sensitive  structures  (solid  curves). 


spatial  importance  distribution,  an  almost  uniform  reduction  of  ~20%  (normalized  to  the 
maximum  sensitive  structure  dose)  in  the  dose  to  the  sensitive  structure  can  be  accomplished. 
If  the  dose  to  the  non-sensitive  structure  normal  tissue  is  not  a  limiting  factor,  the  above  result 
suggests  that  the  target  dose  can  be  escalated  by  ~10%  while  keeping  the  radiation  toxicity  at 
the  current  IMRT  level. 


3.2.  Six-field  IMRT  prostate  treatment 

The  IMRT  plans  obtained  using  three  different  penalty  schemes  are  summarized  in  figures  4 
and  5.  The  isodose  distributions  for  the  three  penalty  schemes  are  shown  in  figure  4.  In  figure  5, 
we  compare  the  DVHs  obtained  using  conventional  inverse  planning  (dotted  curves — obtained 
using  uniform  importance  for  the  prostate  target  and  the  sensitive  structures)  and  the  new 
technique  with  non-uniform  importance  in  both  target  and  sensitive  structures  (solid  curves). 
The  DVHs  of  the  plan  with  non-uniform  importance  to  the  prostate  target  and  uniform 
importance  to  the  sensitive  structure  are  shown  in  dash-dotted  curves.  The  importance  factors 
for  both  target  and  sensitive  structures  were  constructed  based  on  the  computed  capability 
maps  using  the  procedure  described  in  section  2.  It  is  seen  from  figure  4  that  the  dose  sparing 
of  the  rectum  and  bladder  is  significantly  improved  when  the  non-uniform  penalty  scheme 
is  employed.  Remarkably,  the  maximum  dose  to  the  rectum  is  reduced  from  68  to  60  and 
the  fractional  volume  is  dramatically  reduced  in  the  whole  dose  range.  The  reduction  in  the 
low  dose  region  is  more  distinct,  but  the  decrease  in  the  high  dose  region  is  also  evident 
and  perhaps  more  clinically  relevant.  Similarly,  significant  improvements  are  achieved  in  the 
doses  to  bladder  and  femoral  heads.  For  example,  the  fraction  volume  receiving  a  dose  of 
35  is  decreased  from  11%  to  8%  for  bladder.  The  high  dose  tail  in  bladder  is  also  evidently 
suppressed  (from  70  to  66).  Interestingly,  the  dose  uniformity  in  the  prostate  target  is  also 
improved:  the  maximum  dose  of  prostate  target  is  slightly  reduced  from  110  to  107  and  the 
minimum  dose  is  increased  from  85  to  87.  In  optimization,  it  is  generally  true  that  there 
is  no  net  gain  (that  is,  an  improvement  in  the  dose  to  a  structure  is  often  accompanied  by 
dosimetrically  adverse  effect(s)  at  other  points  in  the  same  or  different  structures).  However, 
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Figure  4.  Isodose  distributions  obtained  from  three  different  penalty  schemes:  (a)  uniform 
importance  for  every  structure;  (b)  non-uniform  importance  for  the  prostate  target  and  uniform 
importance  for  other  structures  and  (c)  non-uniform  importance  for  every  structure.  From  the 
centre,  (red,  pink,  yellow  and  blue)  curves  represent  95%,  80%,  50%  and  35%  isodose  curves.  The 
100%  isodose  curve  corresponds  to  a  dose  of  72  Gy  in  this  case. 


one  should  note  that  things  may  be  different  when  different  penalty  schemes  are  used,  or 
more  generally,  when  different  objective  functions  are  used.  The  simultaneous  improvements 
in  both  target  and  sensitive  structures  here  are  a  direct  consequence  of  the  enlarged  solution 
space  when  non-uniform  importance  is  permissible. 

To  better  understand  the  technique,  we  have  also  optimized  the  dose  under  the  condition 
that  the  non-uniform  importance  is  allowed  only  for  the  prostate  target  (i.e.,  the  importance  for 
all  sensitive  structures  is  kept  uniform).  The  corresponding  isodose  distribution  is  shown  in 
figure  4(b)  and  the  DVHs  are  plotted  in  figure  5  as  dotted-dash  curves.  In  this  case,  it  is  found 
that  the  target  dose  homogeneity  is  slightly  improved.  For  example,  the  fraction  receiving 
over  95  is  increased  from  91%  to  96%.  The  maximum  dose  is  reduced  from  110  to  105. 
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Figure  5.  Target  and  sensitive  structure  DVHs  corresponding  to  plans  with  three  different  penalty 
schemes:  (a)  uniform  importance  for  every  structure  (dotted  curves);  (b)  non-uniform  importance 
for  all  structures  (solid  curves)  and  (c)  non-uniform  importance  to  the  prostate  and  uniform 
importance  to  the  sensitive  structure  (dash-dotted  curves). 


By  individualizing  the  importance  for  the  target  voxels,  we  selectively  increased  the  penalties 
to  those  voxels  that  are  likely  to  be  underdosed.  Consequently,  the  target  dose  coverage  is 
improved.  It  is  interesting  to  note  that  there  is  essentially  no  change  in  the  DVHs  of  the 
sensitive  structures.  The  systematic  improvement  in  the  isodose  distributions  when  the  three 
different  penalty  schemes  are  employed  can  also  be  easily  appreciated  from  figure  4. 


%  Volume  %  Volume  %  Volume 
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Figure  6.  The  comparison  of  DVHs  for  paraspinal  tumour  case  between  plan  obtained  from  the 
algorithm  proposed,  denoted  by  solid  line,  and  plan  from  conventional  optimization,  denoted  by 
dotted  dash  line.  The  dose  sparing  of  spinal  cord,  kidney  and  liver  is  evidently  improved.  The 
100%  in  the  x-axis  corresponds  to  a  dose  of  56  Gy  in  this  case. 


33.  Five-field  IMRT  treatment  of  a  paraspinal  tumour 

The  DVHs  obtained  using  the  conventional  and  newly  proposed  IMRT  planning  techniques  are 
plotted  in  figure  6.  The  isodose  distributions  for  the  two  plans  are  shown  in  figure  7.  Similar 
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Figure  7.  Isodose  distributions  for  plans  obtained  with  two  different  penalty  schemes:  (a)  uniform 
importance  for  every  structure  and  (b)  non-uniform  importance  for  every  structure.  From  the 
centre,  (red,  pink,  yellow  and  blue)  curves  represent  95%,  70%,  55%  and  30%  isodose  curves. 


to  the  previous  case,  when  spatially  non-uniform  importance  factors  given  by  equation  (4) 
are  used,  the  target  dose  coverage  and  sensitive  structure  sparing  are  all  improved  in 
comparison  with  the  conventional  IMRT  plan  with  uniform  importance  factors.  For  the  target, 
the  improvement  is  evident  especially  in  the  dose  range  from  90  to  95.  The  fractional  volume 
receiving  a  dose  level  of  90  is  slightly  increased  (from  99.3%  to  99.7%).  A  more  notable  change 
is  found  for  the  fractional  volume  receiving  doses  higher  than  95  (from  95.1%  to  97.5%). 
The  minimum  target  dose  is  increased  from  84  to  90.  The  maximum  target  dose  is,  however, 
slightly  increased  from  105  to  106.5.  By  using  the  a  priori  non-uniform  penalty  scheme,  it 
is  found  that  the  doses  to  the  sensitive  structures  are  dramatically  improved.  As  seen  from 
the  DVHs,  the  spinal  cord  is  better  spared,  especially  in  the  high  dose  region.  For  example, 
the  fractional  volume  receiving  a  dose  above  50  dropped  from  53%  to  24%.  The  fractional 
volumes  of  the  left  kidney  receiving  a  dose  above  30  are  reduced  from  42%  to  19%  and  for 
the  right  kidney,  the  reduction  is  about  7%  (from  17.7%  to  10.7%).  The  improvement  for 
the  liver  is  less  impressive  but  evident.  Once  again,  we  would  like  to  emphasize  that  the  huge 
improvement  in  sensitive  structure  sparing  is  achieved  without  significantly  deteriorating  the 
dose  coverage  of  the  tumour  target. 
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4.  Conclusions 

The  dosimetric  inequivalence  of  the  voxels  is  a  fundamental  feature  of  the  system  and  should 
be  considered  to  obtain  truly  optimal  IMRT  plans.  We  have  introduced  the  concept  of 
dosimetric  capability  to  quantify  the  likelihood  for  a  voxel  to  meet  its  dosimetric  goal  and 
developed  a  new  inverse  planning  formalism  with  the  dosimetric  capability-modulated  voxel- 
dependent  penalty.  Using  the  new  formalism,  we  can  effectively  ‘boost’  those  regions  where 
there  are  potential  problems  in  meeting  the  dosimetric  goals.  With  the  aim  of  obtaining 
a  spatially  more  uniform  target  dose  distribution,  in  this  paper,  we  presented  a  strategy  to 
assign  a  higher  penalty  (or  high  local  importance)  to  those  voxels  with  lower  capabilities. 
However,  it  is  important  to  note  that  other  penalty  schemes  can  also  be  constructed  under 
the  guidance  of  the  capability  map  to  meet  a  different  clinical  requirement.  The  technique 
provides  an  effective  mechanism  to  incorporate  prior  knowledge  of  the  system  into  the  dose 
optimization  process  and  enables  us  to  better  model  the  intra- structural  tradeoff.  Comparison 
with  the  conventional  inverse  planning  technique  indicated  that  the  algorithm  is  capable  of 
generating  much  improved  treatment  plans  with  more  conformal  dose  distributions  that  would 
otherwise  be  unattainable.  Finally,  we  mention  that  the  technique  may  find  applications  in 
dose  optimization  for  many  other  radiation  therapy  modalities,  such  as  prostate  implantation, 
gamma  knife,  micro-MLC  based  stereotactic  radiosurgery  and  other  variants  of  IMRT,  such 
as  tomotherapy  and  intensity-modulated  arc  therapy. 
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